{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "838ab677",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "96d3aee3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def realY(n0,t0,t):\n",
    "    return n0 * np.exp(- t / t0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "085993d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "n0 = 1e16 #cm^-3\n",
    "t0 = 1e-6 #s\n",
    "t = np.arange(0,3 * t0,t0 / 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ed6ba8ca",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x21c71ac4d60>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEQCAYAAACgBo8fAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjmElEQVR4nO3deXhV1b3/8fc3EyFkHiAQEhJmGWQwDAIOpQ6oVaxVq7YOSItUadVeb9V7b6ve3t+911qtbZ1rVaytSCta8VZwxgEVgoxhMswJhCRABiABkqzfHzlSjIEc4CT7nJPP63l4yDl77XO+u8t+slh777XNOYeIiIS+CK8LEBGRwFCgi4iECQW6iEiYUKCLiIQJBbqISJhQoIuIhAlPA93MnjGzMjNb5UfbM83sczOrN7PLm23LMbM3zWyNma02s9w2K1pEJEh5PUJ/DpjkZ9utwA3AX1rY9jzwgHPuFGA0UBaI4kREQomnge6c+wDYfeR7ZtbHzOaZ2RIz+9DMBvrabnbOrQAam7UfBEQ5597ytdvrnNvfTocgIhI0vB6ht+Qp4MfOudOAO4DHWmnfH6g0szlmttTMHjCzyDavUkQkyER5XcCRzCweGAf81cy+fLtTK7tFAWcAI2ialnmJpqmZP7ZNlSIiwSmoAp2mfzFUOueGH8c+xcBS59xGADN7FRiLAl1EOpigmnJxzlUDm8zsCgBrMqyV3RYDKWaW4Xs9EVjdhmWKiAQl83K1RTN7ETgbSAd2AvcA7wKPA92BaGCWc+4/zWwU8AqQAtQBpc65wb7PORd4EDBgCTDNOXewfY9GRMRbnga6iIgETlBNuYiIyInz7KRoenq6y83N9errRURC0pIlSyqccxktbfMs0HNzcykoKPDq60VEQpKZbTnaNk25iIiECQW6iEiYUKCLiIQJBbqISJhQoIuIhIlWA721h1D4bs//nZkVmdkKMxsZ+DJFRKQ1/ozQn+PYD6G4AOjn+zONptv2RUSknbUa6C09hKKZycDzrsmnQLKZdQ9Ugc1t2bWP++YWcqihsfXGIiIdSCDm0LOAbUe8Lva99zVmNs3MCsysoLy8/IS+rKhsL89+vJmXlxSf0P4iIuEqEIFuLbzX4opfzrmnnHP5zrn8jIwW71xt1cSBXRmenczv3y3iQH3DCX2GiEg4CkSgFwPZR7zuCWwPwOe2yMz4l/P6U1JZy+zF21rfQUSkgwhEoL8GXOe72mUsUOWc2xGAzz2qCX3TGZ2byu/fLaLukEbpIiLg32WLLwKfAAPMrNjMpprZdDOb7mvyD2AjUAT8Abi5zar9Z0389Lz+lNUc4IVPj7pOjYhIh9LqaovOuatb2e6AWwJWkZ/G9k5jfN80nliwgWvG5BAXE2yPRxURaV8hfafoT88dQMXeg8xcqFG6iEhIB/ppvVL4xoAMnvxgAzV1h7wuR0TEUyEd6NA0Sq/cf4hnPtrsdSkiIp4K+UAf2jOJ8wZ14+mPNlK1X6N0Eem4Qj7QAW4/tz81dfX84cONXpciIuKZsAj0U7onctGp3Xn2403s3nfQ63JERDwRFoEOcPs5/ag91MCTCzZ4XYqIiCfCJtD7dk3g0uFZzPxkM2U1dV6XIyLS7sIm0AF+8s1+HGpwPPaeRuki0vGEVaDnpnfh8pE9+ctnW9lRVet1OSIi7SqsAh3gx9/si8PxyLtFXpciItKuwi7Qe6bE8d1R2cwu2Ma23fu9LkdEpN2EXaADzPhGP8yM373zhdeliIi0m7AM9MykWL4/phdzlpawsXyv1+WIiLSLsAx0gB+d3YfYqAj+9421XpciItIuwjbQMxI6cfM3+vLm6p0s3FDhdTkiIm0ubAMdYOqEPLKSO/Nfr6+hobHF51aLiISNsA702OhI7rxgIKt3VPPykmKvyxERaVNhHegAF5/anZE5yTzw5jr2Hqj3uhwRkTYT9oFuZvz8W4MorznA4+/rZiMRCV9hH+gAI3JSmDy8B3/4cBPFe3SzkYiEpw4R6AB3ThpIhMH989Z5XYqISJvoMIHeI7kz087ozdzl21myZY/X5YiIBFyHCXSAm87qQ9eETvzy9dU06jJGEQkzHSrQu3SK4l/PH8CybZXMXbHd63JERAKqQwU6wHdG9mRIViL3v7GW2oMNXpcjIhIwHS7QIyKMn180iO1VdTz94UavyxERCZgOF+gAY3qnMWlwJo8v2MDOaj1/VETCQ4cMdIC7LxxIfYPj1/N1GaOIhIcOG+i90rpww/hc/vZ5MatKqrwuR0TkpHXYQAeYMbEvKXEx/PL11TinyxhFJLR16EBPjI3m9nP789mm3cwvLPW6HBGRk9KhAx3g6lHZDMxM4L65q7Uao4iENL8C3cwmmdk6Mysys7ta2J5kZnPNbLmZFZrZlMCX2jaiIiP478uGUlpdx4Nv6gSpiISuVgPdzCKBR4ELgEHA1WY2qFmzW4DVzrlhwNnAg2YWE+Ba28zInBS+P6YXMxduZkVxpdfliIicEH9G6KOBIufcRufcQWAWMLlZGwckmJkB8cBuIKTmL/510gDS4jvxb6+spL6h0etyRESOmz+BngVsO+J1se+9Iz0CnAJsB1YCtzrnvpaKZjbNzArMrKC8vPwES24bibHR3HvxYFaVVDPzky1elyMictz8CXRr4b3m1/idDywDegDDgUfMLPFrOzn3lHMu3zmXn5GRcZyltr0Lh2YycWBXHnxzHSWVtV6XIyJyXPwJ9GIg+4jXPWkaiR9pCjDHNSkCNgEDA1Ni+zEz7rtkMM7BPX9fpWvTRSSk+BPoi4F+ZpbnO9F5FfBaszZbgW8CmFk3YAAQkitfZafGcfu5/Xh7TRnzC3d6XY6IiN9aDXTnXD0wA5gPrAFmO+cKzWy6mU33NfslMM7MVgLvAHc65yraqui2NmV8Hqd0T+Te1wqpqTvkdTkiIn4xr6YV8vPzXUFBgSff7Y+lW/dw2eMLuf70XO69ZLDX5YiIAGBmS5xz+S1t6/B3ih7NiJwUrh3bi5mfbGb5tkqvyxERaZUC/RjuOH8AXRM6cfccXZsuIsFPgX4MX16bvnpHNc8t3Ox1OSIix6RAb8WkIZl8c2BXHnxzPcV79ntdjojIUSnQW2Fm3De56aToPX8v1LXpIhK0FOh+6JkSx0/P7c87a8t4Y5XWTReR4KRA99OU8bkMyUrk56+uomLvAa/LERH5GgW6n6IiI3joyuHUHKjnrpdXaupFRIKOAv049O+WwM/OH8Dba3by14Jir8sREfkKBfpxunF8HmN7p3Lf3EK27dZVLyISPBToxykiwvj1FcOIMONfZi+noVFTLyISHBToJ6BnShz3XjKYRZt38/SHIbmopIiEIQX6CbpsZBaTBmfy4JvrWbOj2utyREQU6CfKzPjvy4aS2Dma219axoH6Bq9LEpEOToF+ElK7xPCry4eytrSG37z1hdfliEgHp0A/SRMHduPq0dk8+cEGFm3a7XU5ItKBKdAD4D8uGkR2Shz/8tdl7D1Q73U5ItJBKdADoEunKB66chgle2r55dzVXpcjIh2UAj1A8nNTuemsPrxUsI23Vuvh0iLS/hToAXT7Of05pXsid89ZwS4t4CUi7UyBHkAxURE8/N3hVNfWc+fLK7SAl4i0KwV6gA3ITODuCwfy9poynvpAd5GKSPtRoLeBG8blcuHQTH41f50uZRSRdqNAbwNmxv3fOZWc1Dhm/OVzyms0ny4ibU+B3kYSYqN5/Psjqa47xK2zlmpVRhFpcwr0NjQwM5H/unQoCzfs4uG313tdjoiEOQV6G7v8tJ58Nz+b379bxHvryrwuR0TCmAK9Hdw3eTCndE/k9peWUbxHTzkSkbahQG8HsdGRPP69kTQ0OG75y1IO1jd6XZKIhCEFejvJTe/CA1cMY/m2Sv77H2u8LkdEwpACvR1NGpLJDybk8dzCzcxdvt3rckQkzCjQ29mdFwzktF4p3PXyCorK9npdjoiEEb8C3cwmmdk6Mysys7uO0uZsM1tmZoVmtiCwZYaP6MgIHr1mJLHRkdz85yXsP6j100UkMFoNdDOLBB4FLgAGAVeb2aBmbZKBx4BLnHODgSsCX2r4yEyK5bdXjeCLsr38+yurtIiXiASEPyP00UCRc26jc+4gMAuY3KzNNcAc59xWAOecLrhuxYR+6fz0nP68srSEJxZoES8ROXn+BHoWsO2I18W+947UH0gxs/fNbImZXdfSB5nZNDMrMLOC8vLyE6s4jMyY2JdLhvXg/nlrmbdqh9fliEiI8yfQrYX3ms8RRAGnARcB5wM/N7P+X9vJuaecc/nOufyMjIzjLjbcmBm/uvxURuYkc9tLy1hZXOV1SSISwvwJ9GIg+4jXPYHm19wVA/Occ/uccxXAB8CwwJQY3mKjI3ny2nzSunTiB88vprSqzuuSRCRE+RPoi4F+ZpZnZjHAVcBrzdr8HTjDzKLMLA4YA+juGT9lJHTimRtGse9AA1NnLtaVLyJyQloNdOdcPTADmE9TSM92zhWa2XQzm+5rswaYB6wAFgFPO+dWtV3Z4WdAZgK/v2YEa3ZUc9usZTRquV0ROU7m1SVz+fn5rqCgwJPvDmbPfbyJe+eu5qazenP3Bad4XY6IBBkzW+Kcy29pW1R7FyPHdv24XDaU7+PJBRvpkx7PlaOyW99JRATd+h90zIx7Lh7EGf3S+bdXVvLJhl1elyQiIUKBHoSiIiN49HsjyUvvwvQXlrCpYp/XJYlICFCgB6nE2Gj+eP0oIiOMG59bTOX+g16XJCJBToEexHLS4njy2tMo2VPL9BeWcKC+weuSRCSIKdCD3KjcVH51+al8unE3t764jPoGPe1IRFqmQA8Bl47I4hffGsS8wlLunrNS16iLSIt02WKIuHFCHlW1h/jtO1+Q2Dma/7joFMxaWmZHRDoqBXoIue2cflTVHuKPH20iqXM0P/lmP69LEpEgokAPIWbGL741iJq6eh56az2JsVHcMD7P67JEJEgo0ENMRIRx/3eGUlN3iHvnriaxczSXjezpdVkiEgR0UjQERUVG8LurRzCuTxr/+rcVvFlY6nVJIhIEFOghKjY6kqeuy2dIVhIzXlzKwg0VXpckIh5ToIew+E5RzJwyity0OH44s4Bl2yq9LklEPKRAD3HJcTH8aeoY0uI7ccOzi1i/s8brkkTEIwr0MNAtMZYXpo4hJjKCa//4mRbzEumgFOhhIictjj9NHUN9g+PKJz/hC43URTocBXoYGZCZwKxpYwG46qlPWbOj2uOKRKQ9KdDDTL9uCcy+6XRioiK4+g+fsrK4yuuSRKSdKNDDUF56F2bfdDrxnaK45g+fsmTLHq9LEpF2oEAPU9mpccy+6XTS4mO47o+f8dlGPcpOJNwp0MNYj+TOzL7pdLond+b6Zxfx0Re6+UgknCnQw1zXxFhmTRtLbloXbpy5mPfWlnldkoi0EQV6B5Ae34kXfziW/t3imfanAq39IhKmFOgdREqXGP78g7EM7pHEzX/+nNdXbPe6JBEJMAV6B5LUOZoXfjCGkTkp/OTFpby4aKvXJYlIACnQO5j4TlE8d+Mozuyfwd1zVvLr+etwTs8oFQkHCvQOKC4miqevy+eqUdk88l4RP529nIP1jV6XJSInSU8s6qCiIiP4n8uGkp0axwPz11FaVccT155GUudor0sTkROkEXoHZmbc8o2+/Oa7wyjYspsrnljI9spar8sSkROkQBe+PaInM6eMZkdlHd9+7GMKt2v9F5FQpEAXAMb1TedvPxpHhBlXPvEJC9aXe12SiBwnBbocNiAzgVduHk9OWhdufG4xsxdv87okETkOfgW6mU0ys3VmVmRmdx2j3SgzazCzywNXorSnzKRYZt80lnF90vjZyyt46K31uqxRJES0GuhmFgk8ClwADAKuNrNBR2l3PzA/0EVK+0qIjeaZG0ZxxWk9+d07X/CTWcvYf7De67JEpBX+jNBHA0XOuY3OuYPALGByC+1+DLwMaPWnMBAdGcGvLj+Vn00awOsrtnPZYwvZskvPKhUJZv4EehZw5GRqse+9w8wsC/g28MSxPsjMpplZgZkVlJfrpFuwMzNuPrsvz00ZzY6qOi7+/Ue8v06/r0WClT+Bbi2813xS9WHgTudcw7E+yDn3lHMu3zmXn5GR4WeJ4rWz+mcwd8YEeiR3Zspzi3n0vSLNq4sEIX8CvRjIPuJ1T6D5Un35wCwz2wxcDjxmZpcGokAJDjlpccy5eRwXn9qDB+avY/oLS6ipO+R1WSJyBH8CfTHQz8zyzCwGuAp47cgGzrk851yucy4X+Btws3Pu1UAXK96Ki4nit1cN5+ffGsTba8q49NGPKSrb63VZIuLTaqA75+qBGTRdvbIGmO2cKzSz6WY2va0LlOBiZkydkMcLU8dQuf8Qlz76sR6YIRIkzKu50Pz8fFdQUODJd0tgbK+sZfoLS1hRXMWPJ/bltnP6ExnR0ikXEQkUM1vinMtvaZvuFJUT9uVDqK/M78nv3y3ihmcXUVZd53VZIh2WAl1OSmx0JPd/51T+57KhLN68m0m//ZB31uz0uiyRDkmBLifNzLh6dA6v/3gC3RJjmTqzgF/8fRV1h455FauIBJgCXQKmb9cEXr1lHFMn5PH8J1u45JGPWFta7XVZIh2GAl0CqlNUJD//1iBm3jia3fsOcckjHzNz4WbdiCTSDhTo0ibO6p/BvNvOYHyfNO55rZCpMwuo2HvA67JEwpoCXdpMenwnnrlhFPdePIiPiiqY9PCHenCGSBtSoEubMjNuGJ/HazPGk9olmuufWcS9rxVqOV6RNqBAl3YxMDOR12ZM4IZxuTy3cDPnP/wBHxdVeF2WSFhRoEu7iY2O5N5LBvPStLFERUTwvac/486/raCqVot8iQSCAl3a3Zjeabxx6xlMP6sPf/u8mHMfWqD1YEQCQIEunoiNjuSuCwby6s3jSYvvxLQ/LeGWv3xOeY2uhBE5UQp08dTQnkm8NmM8d5zXn7cKd3LubxYw5/NiXbcucgIU6OK56MgIZkzsxz9unUDv9C78dPZypjy3mJLKWq9LEwkpCnQJGn27JvDX6eO45+JBfLZxN+c+tIBH3yviQL3WhBHxhwJdgkpkhDFlfB5v3n4mE/qm88D8dZz3mw+0gqOIHxToEpSyU+N46rp8nr9xNJERxtSZBUx5dhEby/XIO5GjUaBLUDuzfwbzbj2Tf7/wFBZv3sP5D3/A/76xlr0HdKepSHMKdAl6MVER/PDM3rx7x1lMHp7FEws28M0H3+fVpSW6GkbkCAp0CRldE2L59RXDmHPzOLolxnLbS8u44olPWFVS5XVpIkFBgS4hZ2ROCq/ePJ77vzOUTRX7uPiRj7j9pWVs3bXf69JEPGVe/ZM1Pz/fFRQUePLdEj6qag/x+PsbePbjTTQ6xzWjc5gxsR8ZCZ28Lk2kTZjZEudcfovbFOgSDkqr6vjdu1/w0uJtdIqKYOqEPH54Zm8SY6O9Lk0koBTo0mFsqtjHg2+u4/UVO0iOi+aWs/ty7em9iI2O9Lo0kYBQoEuHs6qkil/NX8cH68vpnhTL7ef057KRWURF6rSRhLZjBbr+65awNCQriedvHM1ffjiGromx/OzlFZz/8Ae8urSE+oZGr8sTaRMKdAlr4/qk8+rN43ji+6cRFRHBbS8tY+KDC3hx0VatESNhR1Mu0mE0NjreXrOTR98rYnlxFZmJsUw7szdXj86hc4zm2CU0aA5d5AjOOT78ooJH3iti0abdpHWJYeoZeVw7thcJuipGgpwCXeQoFm/ezSPvFrFgfTmJsVHcMC6XKePzSOkS43VpIi1SoIu0YmVxFY+89wXzC3cSFxPJlfnZXD8ul7z0Ll6XJvIVCnQRP63fWcMTCzYwd/l26hsdEwd0Zcr4PMb3TcPMvC5PRIEucrzKaur486db+fNnW6jYe5D+3eKZMj6Pb4/I0k1K4qmTDnQzmwT8FogEnnbO/W+z7d8D7vS93Av8yDm3/FifqUCXUHCgvoG5y3fwzEebWL2jmpS4aK4encN1p+eSmRTrdXnSAZ1UoJtZJLAeOBcoBhYDVzvnVh/RZhywxjm3x8wuAO51zo051ucq0CWUOOdYtGk3z368mTdXlxJhxgVDu3Pt2F6Myk3RdIy0m2MFepQf+48GipxzG30fNguYDBwOdOfcwiPafwr0PPFyRYKPmTGmdxpjeqexbfd+nv9kM7MWb2Pu8u30yejCVaNyuGxkFmnxWuVRvOPPCP1yYJJz7ge+19cCY5xzM47S/g5g4Jftm22bBkwDyMnJOW3Lli0nWb6Id/YfrOf/Vuxg1uJtLNmyh+hI47zBmVw1KpvxfdKJiNCoXQLvZEfoLf1X2eJvATP7BjAVmNDSdufcU8BT0DTl4sd3iwStuJgorsjP5or8bNbvrGHWom3MWVrM/63YQXZqZ77r29YtUXPt0j78GaGfTtOc+Pm+13cDOOf+p1m7U4FXgAucc+tb+2LNoUs4qjvUwPzCUl5avI2FG3YRYTBxYFeuyM/m7AEZdIrSFTJyck52hL4Y6GdmeUAJcBVwTbMvyAHmANf6E+Yi4So2OpLJw7OYPDyLzRX7eKlgG38tKObtNWUkxkZx0andmTw8i9G5qZqSkYDz97LFC4GHabps8Rnn3P8zs+kAzrknzOxp4DvAl5Pi9Uf7DfIljdClozjU0MhHRRX8fWkJ8wt3UnuogR5JsVwyPItvj8hiQGaC1yVKCNGNRSJBYv/Bet5avZNXlpbw4RcVNDQ6BmYmcOmILC4Z1oMeyZ29LlGCnAJdJAhV7D3A/63YwStLS1i2rRIzGNUrlUlDMpk0JFPhLi1SoIsEuc0V+/j7su38Y+UO1u2sAWBYzyQmDenOBUMyydUiYeKjQBcJIRvL9zKvsJR5q0pZUVwFwMDMhMMj9wHdEnRnagemQBcJUSWVtcxbVcr8VaUs3rIb5yAvvQvnDe7GxAFdOa1Xih583cEo0EXCQFlNHW+t3sm8VaV8smEX9Y2OxNgozuyfwcSBXTmrf4aWHugAFOgiYaam7hAffVHBu2vLeG9dORV7D2AGw3omM3FgVyYO7MrgHomamglDCnSRMNbY6CjcXs27a8t4d10ZK4orcQ66JnTirP4ZTOiXzrg+6WQkaPQeDhToIh1Ixd4DLFhXzrvryvhwfTnVdfVA04nV8X3TGd83jdF5acR38udGcQk2CnSRDqqh0bGqpIqPN1TwcVEFizfv4WB9I1ERxvDsZF/ApzM8O5mYKJ1cDQUKdBEBmhYPW7JlDx8VVbCwqIIVJVU4B3ExkYzMSWFUbiqj8lIYkZ1C5xgtJBaMTnZxLhEJE7HRkYdH5QBV+w/xycZdLNxQwaJNu3n4nfU4B9GRxtCsJEblpTI6N5X8XqkkxUV7XL20RiN0ETmsav8hlmzdzaJNe1i0aRcrS6o41OAwgwHdEhiVm8ppvVIYnp1Mr7Q4XUXjAU25iMgJqTvUwNKtlSzevJvFm3fz+ZY97DvYAEBKXDTDs5MZnp3CiJxkhmUnk9RZo/i2pikXETkhsdGRnN4njdP7pAFNJ1nX76xh2bZKlm7dw7Jtlby/vpwvx4V9MrowPDuF4TnJDOuZRP9uCcRGay6+vWiELiInpabuECuKq74S8hV7DwIQFWH065bAkB6JDMlKYkhWIqd0TyQuRmPJE6UpFxFpN845ivfUsqqkilXbq1hVUs2qkip27WsKeTPokxF/OOQHdU9kQGaCli3wk6ZcRKTdmBnZqXFkp8ZxwdDuQFPI76w+wMqSKlaVVFG4vYpPN+7m1WXbD++XHt+JU7onMKBbAgMyEzileyJ9u8ZryuY4KNBFpM2ZGZlJsWQmxXLuoG6H3y+vOcDa0mrWldawtrSGtaXV/OnTLRyobwQgwiA3vQsDMxPo3y2Bvl3j6ds1nrz0LnrgdgsU6CLimYyETmQkZHBGv4zD7zU0Ojbv2sfaHTWsK61mbWkNq0qqeWNV6eGTrxEGOalx9O36z5Dv2zWePhldSIjtuFfaKNBFJKhERhh9MuLpkxHPRad2P/x+7cEGNlbspahsLxvK9lJU3vTzgvVlHGr457nAbomdyE3rQl56F3LTu5Cb1oXeGV3ISY0L++kbBbqIhITOMZEM7pHE4B5JX3m/vqGRrbv3U+QL+Y3l+9hcsY+3Vu88fCIWmk7G9kjqTG563OHAz06NIzsljuzUzmExslegi0hIi4qMoHdGPL0z4jmv2bbqukNsrtjHpop9bK7Yz+ZdTT+/vmIHVbWHvtI2JS6anNQ4eqbGkeML+pzUprDvntQ5JBYvU6CLSNhKjI3m1J7JnNoz+WvbKvcfZNvuWrbu3s+2Pfub/t69n8KSKt4sLP3KNI5Z0/ryPZI70yO5Mz19f2cd8Xdi5yjPl0JQoItIh5QcF0NyXAxDeyZ9bVtDo6O0uo6tu5rCfntlLSV7atleVUthSRVvFe7kYEPjV/aJ7xTVdCVPYmzLfyfFkhoXQ0RE24W+Al1EpJnICCPLN/I+nbSvbW9sdOzad5CSytrDYV9SWUtpVR2l1XV89EUFZTV1NDa7bzM60uiWGMv1p+fywzN7B7xuBbqIyHGKiDDfJZedGJ6d3GKb+oZGKvYepLS6rinoq2oprT5AaVUtXRPb5q5YBbqISBuIiow4PNVCdvt8Z/CfthUREb8o0EVEwoQCXUQkTCjQRUTChAJdRCRMKNBFRMKEAl1EJEwo0EVEwoRnzxQ1s3Jgywnung5UBLAcL+lYglO4HEu4HAfoWL7UyzmX0dIGzwL9ZJhZwdEekhpqdCzBKVyOJVyOA3Qs/tCUi4hImFCgi4iEiVAN9Ke8LiCAdCzBKVyOJVyOA3QsrQrJOXQREfm6UB2hi4hIMwp0EZEwEdSBbmaTzGydmRWZ2V0tbDcz+51v+wozG+lFnf7w41jONrMqM1vm+/MLL+psjZk9Y2ZlZrbqKNtDqU9aO5ZQ6ZNsM3vPzNaYWaGZ3dpCm5DoFz+PJVT6JdbMFpnZct+x3NdCm8D2i3MuKP8AkcAGoDcQAywHBjVrcyHwBmDAWOAzr+s+iWM5G3jd61r9OJYzgZHAqqNsD4k+8fNYQqVPugMjfT8nAOtD+P8r/hxLqPSLAfG+n6OBz4CxbdkvwTxCHw0UOec2OucOArOAyc3aTAaed00+BZLNrHt7F+oHf44lJDjnPgB2H6NJqPSJP8cSEpxzO5xzn/t+rgHWAFnNmoVEv/h5LCHB97/1Xt/LaN+f5lehBLRfgjnQs4BtR7wu5usd60+bYOBvnaf7/nn2hpkNbp/SAi5U+sRfIdUnZpYLjKBpNHikkOuXYxwLhEi/mFmkmS0DyoC3nHNt2i/B/JBoa+G95r/d/GkTDPyp83Oa1mjYa2YXAq8C/dq6sDYQKn3ij5DqEzOLB14GbnPOVTff3MIuQdsvrRxLyPSLc64BGG5mycArZjbEOXfkOZuA9kswj9CL+eqzsnsC20+gTTBotU7nXPWX/zxzzv0DiDaz9PYrMWBCpU9aFUp9YmbRNAXgn51zc1poEjL90tqxhFK/fMk5Vwm8D0xqtimg/RLMgb4Y6GdmeWYWA1wFvNaszWvAdb4zxWOBKufcjvYu1A+tHouZZZqZ+X4eTVPf7Gr3Sk9eqPRJq0KlT3w1/hFY45x76CjNQqJf/DmWEOqXDN/IHDPrDJwDrG3WLKD9ErRTLs65ejObAcyn6SqRZ5xzhWY23bf9CeAfNJ0lLgL2A1O8qvdY/DyWy4EfmVk9UAtc5XynwYOJmb1I01UG6WZWDNxD08mekOoT8OtYQqJPgPHAtcBK33wtwL8BORBy/eLPsYRKv3QHZppZJE2/dGY7515vywzTrf8iImEimKdcRETkOCjQRUTChAJdRCRMKNBFRMKEAl1EJACslcXeTuDzcszsTd9CZat9d84ekwJdRCQwnuPrNw6djOeBB5xzp9C0HlRZazso0EVEAqClxd7MrI+ZzTOzJWb2oZkN9OezzGwQEOWce8v32Xudc/tb20+BLiLSdp4CfuycOw24A3jMz/36A5VmNsfMlprZA74blI4paO8UFREJZb4FxsYBf/WtVADQybftMuA/W9itxDl3Pk3ZfAZNq01uBV4CbqBpWYSjUqCLiLSNCKDSOTe8+QbfomMtLaL2pWJgqXNuI4CZvUrTAzCOGeiachERaQO+ZX83mdkVcPhxc8P83H0xkGJmGb7XE4HVre2kQBcRCQDfYm+fAAPMrNjMpgLfA6aa2XKgED+fVOZbR/0O4B0zW0nTuul/aLUGLc4lIhIeNEIXEQkTCnQRkTChQBcRCRMKdBGRMKFAFxEJEwp0EZEwoUAXEQkT/x/3pXrUsK5bGwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(t,realY(n0,t0,t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e78f2118",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calNewW(wi,s,h):\n",
    "    return wi + h *s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "0ab08cc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(n,t0):\n",
    "    return - n / t0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "9f4f59be",
   "metadata": {},
   "outputs": [],
   "source": [
    "def s_Euler(wi,t0):\n",
    "    return f(wi,t0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "7cc47694",
   "metadata": {},
   "outputs": [],
   "source": [
    "w_Euler = np.zeros(t.size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "a621c34f",
   "metadata": {},
   "outputs": [],
   "source": [
    "w_Euler[0] = n0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "0fe1d2d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "h = t0 / 10.5 # 约在10.5最贴合，步长比这个大，点在线下，步长比这个小，点在线上\n",
    "for i in range(1,t.size):\n",
    "    w_Euler[i] = calNewW(w_Euler[i-1],s_Euler(w_Euler[i-1],t0),h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "d00e95c5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEFCAYAAAACFke6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjnUlEQVR4nO3deXiU5b3/8fc9k0xmgpKwhCWAsgcCooGAWpciWgEVRavFXRFcqra2nupxq4ogYun5edRqKxUUrEKtC3Ut7mItSsK+JoQ9CcgatsySmbl/fyR4YgwQyDLb53VdXBeZLd+n99WvN8/znc9jrLWIiEh8c0S6ABERaXxq9iIiCUDNXkQkAajZi4gkADV7EZEEkBTpAmrTunVr27lz50iXISISUxYsWLDDWptR23NR2ew7d+5Mfn5+pMsQEYkpxpiNh3pOp3FERBKAmr2ISAJQsxcRSQBq9iIiCUDNXkQkATT4NI4xpivwIJBmrb286jEHMB5oDuRba6c39O8FmL2ohMlzCigt85KZ7uGeoVmMzOnQGL9KRCSm1Glnb4yZZozZZoxZXuPxYcaYAmNMkTHmPgBr7Tpr7ZgaH3EJ0AGoAIobovCaZi8q4f63llFS5sUCJWVe7n9rGbMXlTTGrxMRiSl1PY3zMjCs+gPGGCfwHDAcyAauMsZkH+L9WcA8a+3dwC+PrdTDmzynAG8gxBVfJNN+hwHAWxFi8pyCxvh1IiIxpU7N3lo7F9hV4+FBQFHVTj4AzKJyB1+bYmB31d9Dtb3AGHOLMSbfGJO/ffv2upT1A6VlXlrvNZy9NJnHXvYw4j/JOEOVj4uIJLr6XKDtAGyu9nMx0MEY08oY8xcgxxhzf9VzbwFDjTHPAnNr+zBr7RRrba61Njcjo9Zv+x5WZrqHHWmWB8aUs7BHiJ9/5eKR6W5yD3iO+rNEROJNfS7Qmloes9bancBtNR4sB2qex29Q9wzN4v63lrGvWYg/X+Lnm+wgN3yUwu3PG9amrqXzo51xepyNWYKISNSqz86+GOhU7eeOQGn9yjl2I3M68MRlJ9Eh3YMBtg10Yd7vRuaY9mz+w2by++VT9mVZpMoTEYmo+uzs84AexpguQAlwJXB1g1R1jEbmdPjxqOXZ0ObKNhTcXMDiwYvJvC2TlVe7mfx1kUY0RSRh1HX0ciYwD8gyxhQbY8ZYa4PAncAcYBXwurV2ReOVeuxaDGnBwKUD6Xh3R0qmlBK+cC2tFwQ0oikiCcNYayNdw4/k5ubaxoo4HnX7Z4z4B3Tc4WBedpDXzvWzLxU6pHv4+r4hjfI7RUSagjFmgbU2t7bnoirP3hgzAhjRvXv3Rvsd85t7yb8RRsxL5qJ5yfRZn8qr5/mZ31sjmiISv6IqG8da+6619pa0tLRG+x2Z6R5CTph9ZgWP3Ohle3qYX77r5t53U/EV+xrt94qIRFJUNfumcM/QLDzJlSOYJRmWCdf6+MfPKui13kFenzxKp5Riw9F3aktEpD6i6jROUzg4dfN9YFpLDz97MovTmrek4JYCCm8tZNvMbfT8a08+2rdbwWoiEhcS7gLt4Vhr2TptK0X/VUTQF+atMwO83z9AuOrfP55kJ09cdpIavohEpcNdoE240ziHY4yh/Zj2DFo5iFXdwvz802QeesVNx20KVhOR2KZmX4uUzBQmX1TOcxf7aL3XwaPTPVz6VTJJQQWriUhsSrhz9nWV2cJDXm8vKzuXc/WnLi75j4vcgiTevSLSlYmIHL2o2tkbY0YYY6bs2bMn0qV8P7VzwAN/vSjA/1zuw1NhuPUFB2t+s4bQgVqTmkVEolJUNfummLOvq5rBarsGuHC8350Ot2dS8nQJeX3z2PVJzYh/EZHopGmcY1D2VRkFYwrwrvHS7qZ2FFyfyh/mKVhNRCJL0zgNLP2sdHKX5HLCfSewZfpWuGgtbecrWE1Eopea/TFyepx0faIrz/8yTFlqmF+/7eaO2Smk7Tca0RSRqKNpnHrKb+Zl4fUwfH4yl3ydTO+NHmYOCfCfvhrRFJHooZ19PR0MVnvv9AoeHu1lS6swN3+Qwv1vp+LbqGA1EYkOavb1VD1YbUsry8RrfMwaWkGPYgfz+8yn+E/FClYTkYhTs6+nmiOamS08DHsim9NXnUramWkU/aqIRWcv4sDqA5EuVUQSWFSNXla7ecnNa9asiXQ59Wat5btXvqPoN0WEDoQov7klE9tvp3ifRjRFpOEdbvQyqpr9QdE+Z3+0At8F+OK6pbg+3s/GNiGmDQ+wsV1YKZoi0qA0Zx9hrrYuxg3ZyzOX+kg7YHh4hpvLv0wmWK4RTRFpGhq9bCKlZV5KesLqE7xc+ZmLi75xMaAwiZeG+SNdmogkAO3sm0hmugeAcjdMuyDAH0Z5SQrBA695KLyjkOC+YIQrFJF4pmbfRKqPaAKs7Bxmwq0BfFenU/rnUvL65LHzw50RrFBE4plO4zSRH937tmoaZ1hOB/b8ag8FYwpYdsEy2l7XlsLRzfjDtwpWE5GGo2mcKBH2h9k4cSMbJm5kv8sy4zw/eb1CYHTvWxGpG03jxABHioMu47rw7O1htjcPc8c7bn79dgrp+xSsJiL1p2YfZRZ5vEy4zseswX76rncycaqHs5ckUbpbwWoicuyiqtlH020JIyUz3UPYAf86NchDN3nZ1CbMTf9K4aE3U/GuU8MXkWMTVc0+mm5LGCnVp3a2tbA8eZWPVy+ooOtWJ3l989j81GZsKPqus4hIdIuqZi+1B6tdOCGb01cPosW5LVh791oWnrGQ/cv3R7pUEYkhmsaJIdZats3aRtGviwjuCXLigyeyeHgSkz8r1JimiGgaJ14YY2h7VVsGrhxIxhUZbHh0AwcuKsS1yq/734rIYanZxyBXhovsV7OZcX0Yjxd+/4qbUZ+5cFWgMU0RqZW+QRvDPm/vZd5Y+MUXLobnJTNgjZNpw/wUoKkdEfkh7exjWGa6B28KTB8aYNJVXsIG7pvl4fbPUgnuUbCaiPwfNfsYVn1Mc/UJYR4e7eWj04PkLjDMz57Pjnd2RLhCEYkWOo0Tw2qGq2VkeBj0XBYDQsdTcFMByy9ZTsaoDHo804MPSrb/KIRNUzsiiUOjl3EqHAiz6Q+b2Dh+IyGPYfpgH19mVYCpfF7haiLxR6OXCcjhctD5oc7kLsplc1qI0f908ds3Umi5t7Lba2pHJLFEVbNXNk7Da5bdjHGjyvnbuX56bXby+FQP5yxKwtjKWyWKSGKIqmavbJzG0b6lh09ygzx4k5e1mWFu+CiF+15zc1KFJ9KliUgTiapmL43j4NTOjnTLH3/h48UL/HTa4eC3f3Ky6clNhIPhSJcoIo1MzT4B/CBczcD6s5Kp+GdXMi5sxbr71rHw1IXsX6JgNZF4pmmcBLf9ze0U3lFIcGeQTv/diaUXJTP5izUa0RSJQZrGkUPK+HkGg1YOos01bdj0+Ca8F63Bs1zBaiLxRs1eSG6ZTO+Xe/PS6DDJFfDAq26u+cRFSkAjmiLxQt+gle/NbePl25vg8rkuzl2QRM4aJy8NC7Cyi0Y0RWKddvbyvcx0D74U+NvPAky8xkcgCe553c2dn6RSsasi0uWJSD2o2cv3qgerFXUM88hoLx+eEaT/4spgte1vbo9whSJyrNTs5Xs173/bprWH05/tTW7+AFI6pLDi8hUs//ly/Fv8kS5VRI6SRi+lTsLBMMX/r5gNj2zA4Xaw966WPJ6yldI9GtEUiRYavZR6cyQ5OOHeE8hdkouvSxKp47Zx1QuWVmVGI5oiMUDNXo5Kas9UHr28nOnn++m2xcHj0zz8LD8Jn18jmiLRTKOXctRK9nopzoEl3ULcMMfFNZ+mMGhVEi8P90W6NBE5hKja2SviODZkplemZe5qbnnqcj8vXOSj3W4H4172sGHCBsIBBauJRJuoavaKOI4N1Uc0MTCvT4jxt/kJn3c8G36/gQUDF7A3f29kixSRH9BpHDlqNe99e3Aa5/ycDux4ZweFvyxk4akL6fRfnVh2qYvJXypYTSTSNHopDa6irIJ1965jy1+3sK1FmKnD/BScUHlqR/e+FWk8Gr2UJpWcnkzWlCymjg2Dhftnerh+jgu3X8FqIpGi0zjSaP7dqjJY7bKvXJyfn8TJa53MGBpgaTcFq4k0Ne3spdFkpnsIJMOsIQHGX+ujPAV++4abu/6VSmB7INLliSQUNXtpNNWndtZnhnn0Ri/vnV3BySsMedl5fDfzO6LxmpFIPFKzl0ZTM1itXSsPZ/5vNgMX5eLu6mbV1atYfvFyfMX6MpZIY9M0jkSEDVmKnylm/YPrMcmGfb9qxcTUrZTs1YimyLHSNI5EHeM0dPptJwYuH4gvy0Xq49u45gVLxm4Fq4k0BjV7iShPVw+PXnqAacP8nPidgwnTPAz7Ngm/gtVEGpRGLyXiSvd4KTkZlnYNcf3HLq78IoVBq5N4ScFqIg1GO3uJuIPBamXHW5651M9zF/tovdfBo9M9rH94PWG/gtVE6kvNXiKuZrBaXu8Qj93mJzS8ORvHbyQ/J58985SEKlIfOo0jEXe4YLWdH+6k8NZCFp2xiA6/7sCKX7iZPFfBaiJHS6OXEvWC+4Ksu38dpc+VsiMtzLRhflZ2VrCaSE0avZSYlnR8Ej3/1JMpt4SpcMK9f/dw0wcuUn0KVhOpK53GkZgxr4WXvNFw8dfJXPBtMv3WOZlxfoBFPRWsJnIkUbWz120J5XAy0z1UJMGbP63gset97Glm+fXbbu5+PxX/Vn+kyxOJalHV7HVbQjmc6lM7G9uFeex6H7PPqaBvoYO87Dy2ztiqYDWRQ4iqZi9yOLUFqw3+n2wGLcklNTuV1TesZunwpfg26stYIjVpGkfigg1bSp4vYd196wDoOqkr+afD5I8LNaYpCUPTOBL3jMPQ8c6ODFoxiLQz0yj6VRE7RhYSWuvDgsLVJOGp2UtccZ/opt+H/fjHFWHa7TCMf8nDhfOScYY0pimJTaOXEneMMXzQ1ctXY+G6j1O4Yq6LQaudTL0gwGY0pimJSTt7iUuZ6R72NoPnRvp55lIfaQcMj0x3c+M3qYS8oUiXJ9Lk1OwlLlUf01zYM8QDY7180y/ET7805J+ST9lXZZEtUKSJqdlLXKo5ptminYc+L/Wm38f9sAHL4rMXU3hHIcF9QWYvKuGMSZ/R5b73OWPSZ7qIK3FJo5eScEIHQqx/aD3FTxcTbpvE84PLyT+x4vvnFa4msUqjlyLVOJs56f5Ud3K+zmF7uII7Z7m4+T0Xzaqu3WpqR+KRmr0krLTT03jounJm/yTAqauSeOLFVAaucoKF0jJN7Uh8UbOXhNamtYfZZ1Xw6A1edjQPc8c7bn79dgpZDk+kSxNpUGr2ktAOTu0Ut7FMuM7HrMF++q538t/PJlH6YqmC1SRuqNlLQqs+tWMdsGxoMr43OtMitzmFNxey5LwleNfqlI7EPk3jiNTChi1bXtzC2nvWYissXSZ0YcHZhsmfKFhNopemcUSOknEYMm/JZOCKgbQ4twVr/2stuy8phDUKVpPYpGYvchjujm76vtOXWVeGabXbMO5lDyP/rWA1iT0KQhM5AmMMc0708u+xcPWnKYz82kVuQRLThvtZr2A1iRHa2YvUQWa6h/2pMGWEn6cu9+Hxw0OvuBnz71RCBxSsJtFPzV6kDqoHqy3pFuLBMV6+GhDizK8Nef3y2P3Z7ghXKHJ4avYidVAzWK1VWw8n/7U3p3xxCsZhWHLuEgpuLqCirOKInyUSCRq9FKmnkDfEhkc3sPmPm3G1dbHrnlY84S/RiKY0OY1eijQip8dJtye70f/b/niPg+Pu3sLF08Mcd0AjmhI91OxFGkjz3OY8eoOPN88K0H+NkydeTOX0FU68AY1oSuRFVbM3xowwxkzZs2dPpEsROSbF+7y8+5MKHrnRy5aWYW59z81v30jBt8kX6dIkwUVVs7fWvmutvSUtLS3SpYgck8z0yrTM0taWidf4+Nu5fnptdjJxqoeS50uw4ei7RiaJIaqavUisqz6iaR3wSW6Q8bcEcPRvxpo71rB48GLKC8sjXKUkIn2DVqQBHZy6mTyn4PtpnN+NyuK8ZzPZOn0ra3+7lrx+eXQZ14WFQxxM/lTBatI0NHop0oT8W/ysuXMNO97awaZ2YV4c5mdT2zCge99K/Wn0UiRKpLRPoe+bfXntmjBp++CR6W5+PjeZ5KCC1aRx6TSOSAR83NHL12Pgys9cjJjnYkBVsNpaBatJI9HOXiQCMtM9HPDA1AsDTP6FD1cQHnjVzc1zUwnuD0a6PIlDavYiEVB9amdFl8pgtS8Ghjj9G0Ne3zx2fbQrwhVKvFGzF4mAmsFqrdt46P9Cb/p/lYPT42Tp0KWsHr2ail0KVpOGoWkckSgT8oXYOGEjmyZtIrl1MmW/a8XEYKlGNOWINI0jEkOcbiddJ3RlQP4AvC0cNLtnK5e+FKb5fqNgNTlmavYiUer4U45n3LVeXv9pgJPXOpn4ooczlyYpWE2OiUYvRaJY8T4vm0+DBT2D3PRhCmM/TOG0VUlMH6pgNTk62tmLRLGDwWrftbRMutrH9PP9dCt18Pg0D8VPF2ND0XfNTaKTmr1IFPtBsJqBz3OCjL/VjxnUjKLfFLHozEUcWHkgwlVKLNBpHJEoVluw2j2jsjjv6Uy2vbaNNXetIT8nnxMfOpHFQ5OY/JmC1aR2Gr0UiWGBbQGK7ipi26xtFLepDFbb0F7BaolKo5ciccrVxkX2zGxeuS5Ms3J4+BU3oz5PxlWhYDX5IZ3GEYkDn2V6+c9Y+MXnLobPd5GzJomXh/kpULCaVNHOXiQOZKZ78KbA9GEBJl3pxVi4b6aHX36eSnCPgtVEzV4kLlSf2ll9Ypjf3+Tl49OCDMw3zO8znx3v7YhwhRJpavYicaBmsFpGhoeBz/dmwDf9SW6ZzPIRy1l59UoC2wORLlUiRNM4InEuHAizadImNk7YiLO5kx7P9GBerwomf6QxzXijaRyRBOZwOej8cGdyF+Xi6e5h1TWr2HRlAeWbfFhQuFqCULMXSRDN+jSj/9f9ef/CMFkbHEyc6uGcRUkYqzHNRKDRS5EEYpyGN/p6+aKj4cZ/pXDDRymcuiqJl4b5KdWYZlzTzl4kwWSme9iebpk8ysfU4X5O2OZg/EsefrHUQzgYjnR50kjU7EUSzPdjmga+6hfkgTFeVnYNM/xDBwtPW8j+JfsjXaI0AjV7kQRTc0yzWSc3nWdmkf16Nv7NfhbkLmD979cT9muXH080eiki36vYWUHR3UV8N+M7Unun8t29LZm0dZNGNGOERi9FpE6SWyXTe3pvTvrgJPbvDtBs9GZ++kYIV0AjmrFOzV5EfqTV8FaMvy3A5zlBzl+QzIRpHrI3ODSiGcPU7EWkVhu8Xl45P8DEq70EnXDv3z3c9IGLsq0a0YxFavYiUquD978t7BTm4dFe3j0twBnLk5g0LZXtb2+PcHVytNTsRaRW1ZM0K5LgzZ9W8ORNATyZblZctoIVV6zAv9Uf4SqlrvQNWhGpVW33v719VBZD/tyezX/czIZxG9j96W66P9Wdb04KKlgtymn0UkSOyYHVBygYW8Der/eysmuIaef72ZFW2U90/9vI0OiliDS4Zr2akTM3h3cuDtOl2MGEqR7OW6BgtWil0zgicsyMw/B2by9fZhpunOPi2k8qg9WmDVewWrTRzl5E6iUz3cPONMv/XOFnyoV+2u908NhLHq5a5CFcociFaKFmLyL1Uj1Y7T99gzwwtpylPcOc/5GDBQMXsG/hvkiXKKjZi0g91QxWO76Dh26v9qLPW32o+K6CBYMWsPa+tYS8oUiXmtAafBrHGNMVeBBIs9ZeXvXYYGA8sAKYZa394nCfoWkckfhQsbuCtfesZevUrXh6eth2bysmbVewWmOp9zSOMWaaMWabMWZ5jceHGWMKjDFFxpj7AKy166y1Y2p8hAX2A26g+OgPQURiUXKLZHq92It+H/dj//4Kjh9bzJDXQ6T4FazW1Op6GudlYFj1B4wxTuA5YDiQDVxljMk+xPu/stYOB/4bGHdspYpIrGp5XkvG3xpgTm4F5yxK4vGpHvqtdWpEswnVqdlba+cCu2o8PAgoqtrJB4BZwCWHeP/BS/K7gZTaXmOMucUYk2+Myd++XbkbIvFmY7mXmecGePxaHz4X3P2Gm1veS2FvqUY0m0J9LtB2ADZX+7kY6GCMaWWM+QuQY4y5H8AYc5kx5gXgFeBPtX2YtXaKtTbXWpubkZFRj7JEJBodDFZb2yHMIzd6mf2TAINWOXliWirb/r6NaPw2fzypT7M3tTxmrbU7rbW3WWu7WWufqHrwLWvtrdbaUUe6OCsi8al6sFowCWafVcGkMQHcJ7pZeeVKlo9cjr9UwWqNpT7foC0GOlX7uSNQWr9yRCRe1RasdueoLIY8357i/y1mw+83MD97Pt3+2I35/cMKVmtg9Wn2eUAPY0wXoAS4Eri6QaoSkbg0MqdDrU37hN+dQOuRrSm8uZDCmwv57sQQgaF+bIv/m9o5+H45NnUdvZwJzAOyjDHFxpgx1togcCcwB1gFvG6tXdF4pYpIPEvtnsrJn57M25eG6bTFwYRpHobOT8KEFazWEOq0s7fWXnWIxz8APmjQikQkYRmH4Z2eXr5sb7jhIxdXfZ7CoNUKVmsIURWXYIwZYYyZsmfPnkiXIiIRkpnuoex4y9OX+fnzCB8ZZQ7Gvezh2vxUwgEFqx2rqGr21tp3rbW3pKWlRboUEYmQ6sFq32aHeGBsOQuyw5z7qSG/fz57v90b6RJjUlQ1exGRmsFqaZkesmb04qT3TiK0J8TC0xdSdHcRoQMKVjsaui2hiMSM4N4g6+5bR+mfS3F3dbPj3lY8sXuzRjSr6LaEIhIXkpon0fP5npzy5SkcCIY47rYSzp8ZwuNTsNqRqNmLSMxJPzudx8b6ef/UAGcuqwxWy1mjYLXD0T1oRSQmbTrgZeNgmN8rxJgPXdz1lptvewV57TyNaNZGO3sRiUkHg9U2tgsz7nofb54VoP8aJ09MTWXrK1sVrFZDVDV7zdmLSF1VD1YLOeHdn1Qw8eYArh4eVl+/mmUXLsO3yRfhKqNHVDV7zdmLSF3VHNHskO7hrlv6cu7CU+n+THfK5paR1yePkudLsGHt8jV6KSJxybvBS+Ethez+eDfBHDfPnlPO0uT4HtHU6KWIJBxPZw/95vSj/NG2BFZ5ufNpwwXfJLNlV2KOaKrZi0jcMsbweMoW7h/jZUm3EFd86eLhGW4yim3CjWiq2YtIXCst87LnOMufLvXzp5E+Wuw3PDLdzWnvBAn5EidyQc1eROLawRFNgPysEA+M8fJ13yAj5rnIPyWfPV8nxvSfmr2IxLXqI5oABzww85IQ+5/LJOwLs+isRaz51RqC+4IRrLLxRdU3aI0xI4AR3bt3j3QpIhInarv37T1Ds7gopwPB64Osf3A9Jc+WsOOdHWRNyWJuG++PXhsPkzsavRSRhLfnP3soGFNA+epy5vUL8bfBPg5Unf3xJDt54rKTYqLha/RSROQw0n6SxoBFA/j8HMvA5Q4mvphK7urKUz/xEq6mZi8iAjjdTmYMKmfcDT52NQ9z5z/d3Pl2Cmn7DaVlsR+upmYvIlIlM93D5jZhxl/n4++DA/Rb52Tiix4uKvLEfLCamr2ISJWDkzthB3x4agW/H+2lpK3l5286WHr+UrzrY3eHr2YvIlKlZrhaUlc3bd7qSY8/92Dvt3vJ65tH8dPF2FDs7fI1jSMiUge+zT4Kbytk1we7aH5ac4p/l86TazZE1YimpnFEROrJ3cnNSe+dRO+/9aZs9X48ozYy4P0gjlBs3P9WzV5EpI6MMbS9pi2T7giyoGeIy/7t4pHpbjpvcUT9iGZUNXvdqUpEYsGaoJe/XOznfy/zcbzX8PArbkZ9nsyO7dF7ATeqmr3uVCUiseBguNriHiEeGOvly35Bhs938fj0VHZ/sTvC1dUuqpq9iEgsqB6u5k2B6cMCPHVtgBYeF0vOWULBbQUE90RXsFpUBaGJiMSC2sLVRo/KYvAL7Vj/8HqKnypm53s76fmXnvy7gz8qgtU0eiki0sD25u2lYEwBB5YdIK9PiBlDfOxLrXyuMYPVNHopItKEmg9szoD8AXxyXpicVZXBaqetdIKNXLCamr2ISCNwuBy8OsDLwzd62ZYe5rZ33fzmzRRa7I1MsJqavYhII8lM91CaYZlwrY/XhvjpvcnJxKkeLinwYMNNewpdzV5EpJEcnNqxDvhoYJCHbvKyMTPMyNkOFg9ZTPma8iarRc1eRKSR1AxWc3V20+7NLLKmZrF/8X7y++WzafImwsFwo9eiaRwRkQjwl/opvL2Qnf/cyXEDjqP0npY8ub5+wWqaxhERiTIpmSn0fbsv2a9ns3e9l2ZXb+TUd4I4g40TrBZVzV7ZOCKSSIwxtLmiDRNvr2BedpCL57l46G9uTLjhRzSj6hu01tp3gXdzc3NvjnQtIiJNZW2Fl6IL4dveIVrtNdiqbXhDjmhGVbMXEUlEmekeSsq8LOsa+tHjDSWqTuOIiCSi6sFqB3mSndwzNKvBfod29iIiEVZbsFpDB6ap2YuIRIGROR0aNQ1Tp3FERBKAmr2ISAJQsxcRSQBq9iIiCUDNXkQkAURlEJoxZjuwsR4f0RrY0UDlRFK8HAfoWKJVvBxLvBwH1O9YTrTWZtT2RFQ2+/oyxuQfKvktlsTLcYCOJVrFy7HEy3FA4x2LTuOIiCQANXsRkQQQr81+SqQLaCDxchygY4lW8XIs8XIc0EjHEpfn7EVE5IfidWcvIiLVqNmLiCSAmG32xphhxpgCY0yRMea+Wp43xphnqp5faozpH4k666IOxzLYGLPHGLO46s/DkajzSIwx04wx24wxyw/xfCytyZGOJVbWpJMx5nNjzCpjzApjzF21vCYm1qWOxxIr6+I2xsw3xiypOpZxtbymYdfFWhtzfwAnsBboCriAJUB2jddcAHwIGOA04NtI112PYxkMvBfpWutwLGcD/YHlh3g+JtakjscSK2vSHuhf9ffjgcIY/v9KXY4lVtbFAMdV/T0Z+BY4rTHXJVZ39oOAImvtOmttAJgFXFLjNZcAM2ylb4B0Y0z7pi60DupyLDHBWjsX2HWYl8TKmtTlWGKCtXaLtXZh1d/3AauAmqHpMbEudTyWmFD1v/X+qh+Tq/7UnJZp0HWJ1WbfAdhc7edifrzodXlNNKhrnadX/ZPvQ2NMn6YprcHFyprUVUytiTGmM5BD5S6yuphbl8McC8TIuhhjnMaYxcA24GNrbaOuS6zeqcrU8ljN/yrW5TXRoC51LqQy82K/MeYCYDbQo7ELawSxsiZ1EVNrYow5DngT+I21dm/Np2t5S9SuyxGOJWbWxVobAk4xxqQDbxtj+lprq18jatB1idWdfTHQqdrPHYHSY3hNNDhindbavQf/yWet/QBINsa0broSG0ysrMkRxdKaGGOSqWyOr1pr36rlJTGzLkc6llhal4OstWXAF8CwGk816LrEarPPA3oYY7oYY1zAlcA7NV7zDnB91RXt04A91totTV1oHRzxWIwx7Ywxpurvg6hct51NXmn9xcqaHFGsrElVjVOBVdba/3eIl8XEutTlWGJoXTKqdvQYYzzAecDqGi9r0HWJydM41tqgMeZOYA6V0yzTrLUrjDG3VT3/F+ADKq9mFwHlwOhI1Xs4dTyWy4FfGmOCgBe40lZdro8mxpiZVE5DtDbGFAOPUHnhKabWBOp0LDGxJsAZwHXAsqrzwwAPACdAzK1LXY4lVtalPTDdGOOk8j9Ir1tr32vMHqa4BBGRBBCrp3FEROQoqNmLiCQANXsRkQSgZi8ikgDU7EVEGpk5QrDeMXzeCcaYj6pC4VZWfaP4sNTsRUQa38v8+EtT9TEDmGyt7U1lvta2I71BzV5EpJHVFqxnjOlmjPmXMWaBMeYrY0yvunyWMSYbSLLWflz12futteVHep+avYhIZEwBfmWtHQD8Dni+ju/rCZQZY94yxiwyxkyu+nLWYcXkN2hFRGJZVZjbT4B/VKU7AKRUPXcZ8Fgtbyux1g6lsm+fRWXq5ybg78CNVEZJHJKavYhI03MAZdbaU2o+URXwVltg3UHFwCJr7ToAY8xsKm9ucthmr9M4IiJNrCqaeb0x5gr4/haEJ9fx7XlAC2NMRtXPQ4CVR3qTmr2ISCOrCtabB2QZY4qNMWOAa4AxxpglwArqeIe6qhz83wGfGmOWUZl7/9cj1qAgNBGR+KedvYhIAlCzFxFJAGr2IiIJQM1eRCQBqNmLiCQANXsRkQSgZi8ikgD+Px8cdt/RHxMPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(t,realY(n0,t0,t),'m-')\n",
    "plt.scatter(t,w_Euler)\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "c58de2ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "def s_Trape(wi,t0,h):\n",
    "    s1 = f(wi,t0)\n",
    "    wNext = wi + h * s1\n",
    "    s2 = f(wNext,t0)\n",
    "    return(s1 + s2) / 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "f27effb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "w_Trape = np.zeros(t.size)\n",
    "w_Trape[0] = n0\n",
    "h = t0 / 10 # 约在10最贴合，步长比这个大，点在线下方，步长比这个小，点在线上方\n",
    "for i in range(1,t.size):\n",
    "    w_Trape[i] = calNewW(w_Trape[i-1],s_Trape(w_Trape[i-1],t0,h),h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "77960556",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEFCAYAAAACFke6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjbklEQVR4nO3deXRU9d3H8fdvJjOZCUvCLgE0yBIIKltArNaqVQErblVBUQuCik99WlvrU9EKat1pbWurtSgoLrVqiyiiohWQVhEChB0TEtYkYNgSlkwyyczv+SMJjWmAQJbZPq9zOCeZLb/rPXy93PvJ5xprLSIiEt0coV6AiIg0PQ17EZEYoGEvIhIDNOxFRGKAhr2ISAyIC/UC6tK+fXubkpIS6mWIiESUFStW7LHWdqjrubAc9ikpKSxfvjzUyxARiSjGmG1He06ncUREYoCGvYhIDNCwFxGJARr2IiIxQMNeRCQGNHoaxxhzOvAAkGitvbbqMQfwa6A1sNxaO6uxfy7AnMx8ps3PoqDIR3KSl3uHp3LVwC5N8aNERCJKvY7sjTEzjTGFxph1tR4fYYzJMsbkGGPuA7DWbrbWTqj1EVcCXYByIK8xFl7bnMx8Js9eS36RDwvkF/mYPHstczLzm+LHiYhElPqexnkFGFHzAWOME3gOGAmkATcYY9KO8v5UYIm19ufAnSe31GObNj8Lnz/AdYtcdN5jAPCVB5g2P6spfpyISESp17C31i4G9tV6eCiQU3Uk7wf+RuURfF3ygP1VXwfqeoEx5nZjzHJjzPLdu3fXZ1nfUlDko/0Bw/dWu3jkFS+jvnThDFQ+LiIS6xpygbYLsKPG93lAF2NMO2PMC8BAY8zkqudmA8ONMX8EFtf1Ydba6dbadGtteocOdf627zElJ3nZk2iZPLGElb0C/PBfbqbO8pB+2HvCnyUiEm0acoHW1PGYtdbuBSbVerAEqH0ev1HdOzyVybPXcrBFgD9fWcZXaRX86JN4/ud5Q25CLikPpeD0OptyCSIiYashR/Z5QLca33cFChq2nJN31cAuPHHNmXRJ8mKAwiFuzLweJN/amR1P72B5/+UULS4K1fJEREKqIUf2GUAvY0x3IB8YA9zYKKs6SVcN7PLfUcvzoeOYjmTdnsWq760i+c5kNtzgYdoXOYpoikjMqG/08k1gCZBqjMkzxkyw1lYAdwHzgY3A29ba9U231JPX5vttGLJmCF1/3pX8vxQQvDyXdiv8imiKSMww1tpQr+G/pKen26aqOB79Pwu44m3ostfBl2kV/PX7ZRxKgC5JXr6476Im+ZkiIs3BGLPCWpte13Nh1WdvjBkFjOrZs2eT/YxlrX2sGAeXL3Fx+VcuztiawOsXl5HRRxFNEYleYdWNY62da629PTExscl+RnKSl4o4mPPdch4a52NP6yD/876H/5ubQFlBWZP9XBGRUAqrYd8c7h2eitdVGcHM62B59OZS/n5xOX22OFiWtoyClwoIx1NbIiINEVancZpDdeqmujCtc1svFz+dytmt2pB9WzbZt2VT+GYhqS+mMr94n4rVRCQqxNwF2mOxQcvOl3aSe28uFf4gf/+unw8H+LFV//7xupw8cc2ZGvgiEpaOdYE25k7jHItxGJJvT2bI+iFkpQS57lMXv3rdQ5fdKlYTkcimYV8HT1cPT11Rwp9HldKhyMHDr3i58t8qVhORyBVz5+zrK7mNl6VpPjaklHDjZ/Fc/YWb9Kw45l4f6pWJiJy4sDqyN8aMMsZMLy4uDvVSjqR2DibAX0aV8fsfltKyzHDnCw5yfpFDoKTOpmYRkbAUVsO+OXL29VW7WG13uhvzQQ+Sb+tM3m/zyDgrg/2L9h/3c0REwoHSOCdh/8L9ZN2WRWluKZ1v78zGsV6mfaliNREJLaVxGlmbCyuL1br9ohsFL+3EXp5Lh+UqVhOR8KVhf5KcCU56TOvBC5OCHIq33P0PD3e8H0+rEkU0RST8KI3TQMta+Vg+rrJYbdSSymK1Ny4uY2lfRTRFJHzoyL6BkpO8BJzw3nnlTBnnozApyKS5Hn75fgKleaWhXp6ICKBh32A1i9UKOlgevamUdy4pJ3Wrg4y0DAr+UoANht9FcBGJLRr2DVQ7opnc1sslT6UxbP1QWg1pRfakbFZ/fzUlOSWhXqqIxLCwil7WuHnJbZs2bQr1chrMWsuumbvIuScHW2Y5PKktj7cvJO+gIpoi0viOFb0Mq2FfLdxz9ieqrKCMxWPX4Fp0mM2nBJg5soy8jlYtmiLSqJSzD7H45HgeuvQgz11RSvsDDh6a5eWqf7ko9ymiKSLNQ9HLZlJQ7CO/L1XFam6u+tLNkKw4Zo7UrRBFpOnpyL6ZJCd5ATjshRcv9/Pba0vx+OGBNzzk/DyHwGEVq4lI09GwbyY1I5oAa3sEeHSSn/Jrk8j7XR4ZZ2aw/zMVq4lI09BpnGZS+9631Wmc4QO7ULS4iKyJWay+eDWdJ3bm65u9PK1iNRFpRErjhImAL8DWh7ey/Tc7KPZaZl1aRmavylM7Su2ISH0ojRMBnF4nPZ7swfN3BjmQEOSnsz3c+V48rQ6rWE1EGk6nccLM8hY+Vt4Cly11ccWXLvpVFat9laZiNRE5eWF1ZB9OtyUMlepitbnfqSxW29k2yB0feLhvTgKl21WsJiInJ6yGfTjdljBUaqZ2dra3PD62lLeGl9N7h4OMfhnkP5+vYjUROWFhNeyl7mK14U+kMWzDUFqf05pNP97EqgtWUZKtYjURqT+lcSKItZZds3aR+7NcAr4A3R/uzsqLHEz7LFsxTRE5ZhpHF2gjiDGGzuM603Z4Wzb9eBOb79vMgVOCOEeUYTv95/63gAa+iHyLTuNEoPjO8Zwx+wz+OjZI4kGYOsvDDxe7cFUopikiddORfQT7tKuPLybADQvcjFriZnBVsVouimmKyLfpyD6CJSd5OeyFl37g5zfXleKugPvf8HDb4gQqDlWEenkiEkY07CNYzZjmutMD/OpWH5+nBzjnK0PGGRns+2RfiFcoIuFCp3EiWO1ytXadvAycnsrAwy3ImpjFmuFrOGXcKfR4pgfzthb+VwmbLuKKxA5FL6NUoDTAtke2sf3p7QTbOHnxAh9LepQfeV7laiLRR0VoMcjpcXL646czOGMwu1wV3PF3Nz9+N57EQwZQakck1oTVsFc3TuNrNbAVU8aW8M75fgbkOnlshpdz18aBhYIipXZEYkVYDXt14zSNTu28zDunnCnjfRS0C3Lbh/Hc83Y8adYb6qWJSDMJq2EvTaM6tbOzneWJsaW8dkkZvQqc/OKPTvL+mKdiNZEYoGEfA2qWq2Hg64tclM3uTtvzk8j5SQ6Z52dy+OvDoV6miDQhpXFimLWWb177hpy7cwgcDpAyNYXMi51MW6BiNZFIpDSO1MkYwym3nMLQjUNpf2V7tjywhYNXZBP3dRmW/xSrzcnMD/VSRaSBNOwFdyc3/d7uxxs3BWl1CKa86uHaz124yhXRFIkW+g1aOeKfXXx8MRHGLHBz+VduBmdXFqvlqFhNJOLpyF6OSE7yUuKBmZf5eXq0j7gAPPCGlzs+T6DioIrVRCKZhr0cUbNYbUNKkF/d6mPB0ArOXlpZrLb3470hXqGInCydxpEjaherte/oZdALqQwqbUnWhCzWjlxLp1s60fOZnszbrmI1kUii6KXUS7AsyLZHt7H9ye0EWjt46UIfX/Yoh8qqHRWriYQBRS+lwRzxDrr/ujuDlw9mp7eC2//h5q45KlYTiRQa9nJCWvZvyZQbS3jrAj9nbXby+EtezlujYjWRcKdz9nLCTmnr5aOzfazsVcH4j+OZ+FE8wzbE8fH1JtRLE5GjCKsje1UcR4bq1M43bS1P3VDKrEvL6LHTwT1/dJL3hzxsIPyuA4nEurAa9qo4jgy1i9WyL3Thn9Oddhe2IefuHDLPy+TwBhWriYQTpXGk0VhrKfxrIZt+uonAwQCn/eo0Vg2PU7GaSDNRGkeahTGGTmM7MXTDUDpc04GtU7ZyaFQ2ro0qVhMJNQ17aXTujm7S3kzjtZuDtCiBKa95GL3QhVvFaiIhozSONJkFyT6WTIDRC92MXOZm4KY4Xh5RRraK1USanY7spclUF6u9PNLPU2N8OCxMftPLpEUJVBxQsZpIc9KwlyZTs1ht42mVxWr/PLuCoRmGjH4Z7J2nYjWR5qJhL02mZkTTAB06eEn/c18GfzkIZ6KTtZevZcNNG/Dv8Yd6qSJRT9FLCYmgP8i2x7ex/fHtxCXGUXRPOx63BRQUK6IpcrIUvZSw43A76P5QdwavGIyvk4MWk3fxw5lBEg8aRTRFmoCGvYRUyzNb8tCNPt68sIx+WyuL1b63Kg6fXxFNkcak6KWEXP4BH3lDYWWvAOM/jmf8/HjO3hjHKyNKQ700kaihI3sJueQkLwC721ieHlPKy8PLSNnl4NGXvex4ZoeK1UQagYa9hFzNiCYGPh9QwSOTyjDntCT3nlxWfmclh9YdCu0iRSKcTuNIyNW+921ykpd7R6dy8YBkCt8qJOd/c1gxaAWnPXAaq0aqWE3kZCh6KWHPv8dPzk9zKPxrIfkdgrw0oowtyUFA974VqUnRS4lo7vZu0t5I49Vbgnh98ODrHsYscKtYTeQE6DSORIyFnX18ORGuX+RmRIaLQZuczBxRRpaK1USOK6yO7HVbQjmW5CQvpfHw6nA/T9zgI2jgvr95uXNhAhXFKlYTOZawGva6LaEcS83UTtapQaaM9/HJORUMWW5YlraMPXP3hHiFIuErrIa9yLHUVaw29Lm+DF46CFc7F+uuWMeGGzfg361iNZHalMaRqBD0B9n+1Ha2/XobztZOej3biyV9ypn2iWKaEjuUxpGo53A7SHkwhfTMdLw9vWwcu5Edo7PwbS/V/W9F0LCXKNOiXwsGfTGIeT8I0nubg8dmeLlgVRzGKqYpsU3RS4k6xmn4+xk+FnU1jPs4nnHz4xm2ofL+twWKaUqM0pG9RKXkJC+7kyzTRpcyc0QZpxY6+PXLXq5f4yVYEQz18kSanYa9RKUjMU0Di/tXcP8EHxtODzLyIweZ52RyaI2K1SS2aNhLVKod02zRzUPKm6mkvZVG6bZSVgxewZapWwiWBZmTmc+5Ty6g+33zOPfJBbqIK1FJ0UuJOeV7y8n5WQ7fvPYNgdPd/Pb8g2zo9J/fwFW5mkQqRS9FanC1c9H31b6c+eGZHNxbxi9ecXPDZ27cVb+LpdSORCMNe4lZ7Ua2Y/L4EhYOrGD4chePzvTSd2vlX4mCIqV2JLpo2EtMa9vJy2uX+nn8Rh9BB/zyLS/jP3Jzerw31EsTaVQa9hLTqlM72d2CPDjex7yz/Zy3No4Hnnexe87uUC9PpNFo2EtMq5naqXDBl1fFUfLaqbTq4mH91etZf/16/N+oWE0in9I4InUIlgfZ8fQOtj6yFWdLJz1/35Ov+lWoWE3CmtI4IifI4XJw2gOnkb4qnYTUBL6+5Wvyr1exmkQuDXuRY2jRtwUD/zWQuaOC9NheWaz2/ZUqVpPIoyI0keMwTsPsNB+fd6ksVrv503iGblSxmkQWHdmL1ENykpc9iZbfXF/KS5eV0XV3ZbHamFVeguUqVpPwp2EvUg81i9X+fWYF90/0sa5XkOHzHaw8eyUHMw+Geokix6RhL1IPtYvVWnb10P2vfej3936UFZSxYsgKNj+wmUBpQMVqEpYUvRRpoPJ95eTek8uuV3YRSHHxzPcOsf4UFatJ81P0UqQJudq66PNyH86afxYHivzcM8vN2H+6iVexmoQRDXuRRtL20rZMHl/CZ4Mr+P6KOB6b4aXfFiegYjUJvbAa9saYUcaY6cXFxaFeishJadfRyxsX+3l8bCn+OLj3bQ8T5rnp4VaxmoRWWA17a+1ca+3tiYmJoV6KyEmpTu3kdA0ydbyP98/x8531cdz/vIvd/1CxmoROWA17kUj3rWK1OFh6RRwlfz2N1qd6WX/tetZdu46yXWWhXqbEIKVxRJpBsCJI3m/z2DJ1C84EJz2e6cHSswIqVpNGpTSOSIg54hyc+stTGbJ6CC36tSBrfBYF12ZRtk3FatI8NOxFmlFCagIDPh/A+1cE6Z7v4NEZXi5eoWI1aXoqQhNpZsZheLevj0VdDOM+dnPTP+M5e2McM1WsJk1IR/YiIZCc5GVfa8sz15Ux/QdldN7r4JFXvNyQqWI1aRoa9iIhULNY7cszKrh/Yglrege59BMHK4eu5OBKFatJ49KwFwmB2sVqrbp46fFGH/rN7od/l58VQ1ewefJmAr5AqJcqUULRS5EwU76/nNx7c9k1Yxfe3l4Kf9mOJwu3K6Ipx6XopUgEcbVx0eelPpz16VkcOlROqwl5XPhOgPgyRTTl5GnYi4Spthe35dd3+JmfXs5FKyuL1c7c7FREU06Khr1IGNtW4uPN7/t57KZSSt1wzzseJs5zU7xTEU05MRr2ImEsOamyLTO3S5Cp43y89x0/wzbE8eTMBArfKSQcr7lJeNKwFwljRyKaQEUcvPvdcp6c4MfTzcOG6zew/pr1lO1UsZocn36DViSMVadups3POpLGuWt0Khc935m83+WxdcpWlvVdRs9nerJ0gIrV5OgUvRSJYCWbSsiamEXx4mI2dg8w49Iy9iRV/p3WvW9jj6KXIlEqoVcCAxYO4L2rgqTkO3hsppdLlsdhgipWk2/TaRyRCGcchjmpPj7vbPjRfDdjP6sqVhupYjX5Dx3Zi0SB6mK1311bxl8uL6XTfgcPv+Jl7AovQb+K1UTDXiQq1CxWW9IvwP0TSliVGuTifzpYMWQFB5YfCPUSJcQ07EWiQO1itdZdvPR6vQ9nvHcG5XvKWXn2SnL/L5dAiYrVYpXSOCJRrqK4gtx7c9n54k68Pb2kvpTKotaHvxXnVEwzOiiNIxLD4hLjSJ2eSv/P+mODllUXrGLNxI3s/can+9/GEA17kRjR5qI2DFk7hH+fZzkv08ljM7z0z6387VzFNKOfhr1IDHEmOJlxbgmP3lSKLx5+9ncPd8yNp1UJFBQpphnNlLMXiTHJSV4242PqOB+XL3ExaomLflsS+PBKi7UWY0yolyhNQEf2IjGmOqYZcMJ755UzdZyPvW0sY/7mYN2V6yjLV7FaNNKwF4kxtWOa9PLQ5r3e9PhtD/b/cz/L0pZR8GIB1lrmZOZz7pML6H7fPM59coEu4kYwRS9F5Ahfro+s27IoWlhEebqXR4YWsaNVxZHnVa4W3hS9FJF68fbw0v+z/vR+sTeBtSX8arqbEcsqi9VAqZ1IpmEvIt9ijCF5YjKTJ/jYkBJgzMJ4HnzdQ5fdlRduldqJTBr2IlKnhG4e/nBNGX8eVUr74spitav+7aJbS2+olyYnQcNeROp07/BUvG4nS9Mqi9WW9Qlw1Rdupr7i4cAyFatFGg17EalTzdTO4QSYd7ODQ39IJqHMsPKcleTck6NitQiiNI6InJCKAxVs/uVmCl4owHO6h9SXUvk8qUTFamFAaRwRaTRxrePo/efeDFg0AOMwrL5oNWsnbGSfitXCmoa9iJyUpO8lkb4mncXnW85d5eTxl7wM2KRitXClYS8iJ83pdfLyOSU8cnMph7xw92wPk96Pp9VhRTTDjYrQRKRBkpO8bMXHQz/ycdlSF1d86aLf1gQ+GqVitXCiI3sRaZCaxWpzv1NZrLa7reX6tx2sHbWW0h2loV6ioGEvIg1Uu1jN9PTQdk5vev6+J0ULi8jol0H+C/nYYPgl/2JJo0cvjTGnAw8Aidbaa6seuwD4NbAe+Ju1dtGxPkPRS5Ho4NvsI+v2LIo+KyLx/ER2/DyJpzZuUUSziTQ4emmMmWmMKTTGrKv1+AhjTJYxJscYcx+AtXaztXZCrY+wwCHAA+Sd+CaISCTynu6l/6f9SZ2RStHKg3iu3cpZ8yswQUU0m1t9T+O8Aoyo+YAxxgk8B4wE0oAbjDFpR3n/v6y1I4FfAg+f3FJFJBIZY+h8a2ee+t8K1nUPMHqRmwdf89Ct0KGIZjOq17C31i4G9tV6eCiQU3Uk7wf+Blx5lPdXFaSyH4iv6zXGmNuNMcuNMct3795dr8WLSOTICvp49uoynruilLYHHEyd5eGaxS4K9yii2RwacoG2C7Cjxvd5QBdjTDtjzAvAQGPMZABjzDXGmL8ArwF/quvDrLXTrbXp1tr0Dh06NGBZIhKOkpO8YCCjb4D7J5awtG8FVyxx89hrCRQvKQ718qJeQ4Z9XeFZa63da62dZK3tYa19ourB2dbaO6y1o493cVZEolN1RBPgsBdevNzPn8b4ae9wkXluJpvu3kTgsIrVmkpDfqkqD+hW4/uuQEHDliMi0ao6dVOzMO2m0amcP70TmydvJv8P+ex9by+9X+zN4nY+Fas1soYM+wyglzGmO5APjAFubJRViUhUumpglzqHdu8/9abj6I5kTchizSVrWN+/gv0XlGE9/0ntVL9fTk59o5dvAkuAVGNMnjFmgrW2ArgLmA9sBN621q5vuqWKSDRL+m4S6avT+fx7lmFrnDw2w8ugbBWrNZZ6Hdlba284yuMfAh826opEJGY5vU5eGVbCghQHt37k5ifveliWWsHrl5RRgFI7DRFWRWjGmFHAqJ49e4Z6KSISIslJXrbh45FbShm5zMWVX7hI26ZitYYKq24ca+1ca+3tiYmJoV6KiIRIzWK1D84pZ8p4H9+0t1z3joO1l62ldLuK1U5GWA17EZHaxWqOHh7azelNz2d7UvSvqmK151SsdqJ0D1oRiRi+rT6yb89m/6f7STyvqlgtS8Vq1XQPWhGJCt4UL2fNP4vUl1PZv/og3uu2MuBjFavVh4a9iEQUYwydx1UWq63uEeC6z91MedXDqd+oWO1YNOxFJCJlB3z86eoy/nRVKW0OGabO8vBDFasdVVhFL0VE6is5yUt+kY/lqQE2nupjzAI3o5a4GZbjovjqYhLPVaqvprA6sjfGjDLGTC8uVgOeiBxb7WK1GT/w8+wNftrFucj8biabfrKJikMVIV5l+FAaR0Qi1pzM/P8qTLu8Vye23L+F/D/lE98tntTpqSzuGBvFasdK42jYi0hUKv6imK8nfI0vy8eSswK8fkEph72Vz3ldTp645syoG/iKXopIzEk8N5H0VeksusAydK2Dx2d4Sc+K3WI1DXsRiVpOj5NZZ5fw8I9K2d/SctccD3e9G0/iIUNBUWyldpTGEZGolpzkZXtVsdqIZS6u/reLvtu8fHx5bBWr6cheRKJadWon6IAPh5Xz4HgfOztafvgPB2uGr8G3NTaO8DXsRSSq1S5Wc/bw0P7d3vR6rhcHlhwg44wM8p7NwwbCL6zSmMIqjVOjz/62TZs2hXo5IhLlSreXkn1HNvs+3kfrc1qTOiOVT0uLIjamqeiliMhRWGv55vVvyLk7h/KDAd4/18/76X4ClcGdiIppKnopInIUxhhOufkUhm4Yyto+Qa5c5GLqLA+n7aocj9ES09SwFxEB3J3c/O6yEp69upRWPsOUVz1ct8iFq5yoiGkqeikiUiU5ycvK3j6+PtXH6IVufrDUzeDsOOZeF+qVNZyO7EVEqlTHNEs88PJIP0+N9hFn4bbpDrJ/nE3FgcgtVtORvYhIleqLsNVpnAMD4on7SS+6vlNG3h/y2Dt3L73/0pt2I9vVWcIWzhdxlcYREamH4q+KyZqQRcmGEvw/aMX9qbvZ4woceT4cUjtK44iINFDisETSV6Zz2oOn4fzoAFNeiGfIRidUHS+He2pHw15EpJ4c8Q66P9Kdh37kY0/rID9+38NP3o0n6WBlv044p3Y07EVETlCwt4dHby7lrQv8nLHFyeMzvJy/Oo7kRG+ol3ZUYTXsdVtCEYkE9w5PJT7eyUdnl/OrW31s7xjk1o/jeejdFvg2h+fRfVgNe2vtXGvt7YmJulGwiISvmuVqu9tYXp9kKLm/I56v/WScmcGO3+8Iu2I1pXFERBpJaV4p2ZOy2TdvH63ObkWfGX341N98xWpK44iINANPVw9nzj2Tvm/0xZfjY9mADBb/bAO79vqwQH6Rj8mz1zInM7/Z16ZhLyLSiIwxdLqxE0M3DmVNmuWKz108NMtL952hLVbTsBcRaQLuDm5+P7KE319TSksfPPiah+sXunCHqFhNdQkiIk0kOcnLql4+sk71cf1CN5ctczNoUxwfXNv8a9GRvYhIE6kuVvPFw6wRfp4c48MBTHzRQdakLCqKm69YTcNeRKSJ1L7/7cH+8bg/6EXXe7qy88WdLOu3jL3z9jbLWhS9FBEJgQPLDpA1IYvD6w7T8caObJrQgqeX5TYooqnopYhImGk9tDWDVwwm5aEUvnm7EMeozXT90o+1TRPR1LAXEQkRh9tBytQUnv1xkMLEIHfOrSxWM8HGj2iGVRrHGDMKGNWzZ89QL0VEpNms8vhYdRNcsiKONocMtuowvDEjmmF1ZK9uHBGJRclJXqwDPhlSwVsXln/r8cYSVsNeRCQWVUc0a/K6nNw7PLXRfkZYncYREYlFte992xSFaRr2IiJh4KqBXZr0/rU6jSMiEgM07EVEYoCGvYhIDNCwFxGJARr2IiIxICyL0Iwxu4FtDfiI9sCeRlpOKEXLdoC2JVxFy7ZEy3ZAw7blNGtth7qeCMth31DGmOVHa36LJNGyHaBtCVfRsi3Rsh3QdNui0zgiIjFAw15EJAZE67CfHuoFNJJo2Q7QtoSraNmWaNkOaKJticpz9iIi8m3RemQvIiI1aNiLiMSAiB32xpgRxpgsY0yOMea+Op43xphnq55fY4wZFIp11kc9tuUCY0yxMWZV1Z8poVjn8RhjZhpjCo0x647yfCTtk+NtS6Tsk27GmIXGmI3GmPXGmJ/W8ZqI2C/13JZI2S8eY8wyY8zqqm15uI7XNO5+sdZG3B/ACeQCpwNuYDWQVus1lwEfAQYYBiwN9bobsC0XAB+Eeq312JbzgUHAuqM8HxH7pJ7bEin7pDMwqOrrVkB2BP9dqc+2RMp+MUDLqq9dwFJgWFPul0g9sh8K5FhrN1tr/cDfgCtrveZK4FVb6SsgyRjTubkXWg/12ZaIYK1dDOw7xksiZZ/UZ1sigrV2p7V2ZdXXB4GNQO3S9IjYL/XclohQ9d/6UNW3rqo/tdMyjbpfInXYdwF21Pg+j//e6fV5TTio7zrPqfon30fGmH7Ns7RGFyn7pL4iap8YY1KAgVQeRdYUcfvlGNsCEbJfjDFOY8wqoBD41FrbpPslUu9UZep4rPb/FevzmnBQn3WupLLz4pAx5jJgDtCrqRfWBCJln9RHRO0TY0xL4B/A3dbaA7WfruMtYbtfjrMtEbNfrLUBYIAxJgl41xhzhrW25jWiRt0vkXpknwd0q/F9V6DgJF4TDo67Tmvtgep/8llrPwRcxpj2zbfERhMp++S4ImmfGGNcVA7HN6y1s+t4ScTsl+NtSyTtl2rW2iJgETCi1lONul8iddhnAL2MMd2NMW5gDPB+rde8D9xSdUV7GFBsrd3Z3Auth+NuizHmFGOMqfp6KJX7bW+zr7ThImWfHFek7JOqNc4ANlprnznKyyJiv9RnWyJov3SoOqLHGOMFLga+rvWyRt0vEXkax1pbYYy5C5hPZZplprV2vTFmUtXzLwAfUnk1OwcoAcaHar3HUs9tuRa40xhTAfiAMbbqcn04Mca8SWUaor0xJg+YSuWFp4jaJ1CvbYmIfQKcC9wMrK06PwxwP3AqRNx+qc+2RMp+6QzMMsY4qfwf0tvW2g+acoapLkFEJAZE6mkcERE5ARr2IiIxQMNeRCQGaNiLiMQADXsRkSZmjlOsdxKfd6ox5pOqUrgNVb9RfEwa9iIiTe8V/vuXphriVWCatbYvlf1ahcd7g4a9iEgTq6tYzxjTwxjzsTFmhTHmX8aYPvX5LGNMGhBnrf206rMPWWtLjvc+DXsRkdCYDvyvtXYw8Avg+Xq+rzdQZIyZbYzJNMZMq/rlrGOKyN+gFRGJZFVlbt8B3qlqdwCIr3ruGuCROt6Wb60dTuXc/i6VrZ/bgbeAcVRWSRyVhr2ISPNzAEXW2gG1n6gqeKursK5aHpBprd0MYIyZQ+XNTY457HUaR0SkmVVVM28xxlwHR25B2L+eb88A2hhjOlR9fxGw4Xhv0rAXEWliVcV6S4BUY0yeMWYCMBaYYIxZDaynnneoq+rB/wXwmTFmLZW99y8edw0qQhMRiX46shcRiQEa9iIiMUDDXkQkBmjYi4jEAA17EZEYoGEvIhIDNOxFRGLA/wNbim5UxQhegAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(t,realY(n0,t0,t),'m-')\n",
    "plt.scatter(t,w_Trape)\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "id": "5a3aa66a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def RK4(wi,t0,h):\n",
    "    s1 = f(wi,t0)\n",
    "    s2 = f(wi + s1 * h / 2,t0)\n",
    "    s3 = f(wi + s2 * h / 2,t0)\n",
    "    s4 = f(wi + s3 * h / 2,t0)\n",
    "    return(s1 +  2 * s2 + 2 * s3 + s4) / 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "id": "36b40d35",
   "metadata": {},
   "outputs": [],
   "source": [
    "w_RK4 = np.zeros(t.size)\n",
    "w_RK4[0] = n0\n",
    "h = t0 / 10.1 # 约在10.1最贴合，步长比这个大，点在线下方，步长比这个小，点在线上方\n",
    "for i in range(1,t.size):\n",
    "    w_RK4[i] = calNewW(w_RK4[i-1],RK4(w_RK4[i-1],t0,h),h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "id": "ac97dbbc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEFCAYAAAACFke6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjZUlEQVR4nO3deXhU9d3+8fd3JjOZCUvCLgEUZAkElS0gVmvRqoAVt6rgWhRUfPRpba2/ilZQ685T29raWhQUl1qXIoqoaAWkVYQAYceEhDUJELaEJZNMMvP9/ZGBxjRAIMts9+u6uK5kluR7PBcfD+fccx9jrUVERGKbI9wLEBGRxqdhLyISBzTsRUTigIa9iEgc0LAXEYkDCeFeQG3atm1ru3btGu5liIhElWXLlu221rar7bmIHPZdu3Zl6dKl4V6GiEhUMcZsOdpzOo0jIhIHNOxFROKAhr2ISBzQsBcRiQMa9iIicaDB0zjGmNOBh4Bka+01occcwG+AlsBSa+2Mhv69ALOyCpgyN5vCYh+pKV7uH57GlQM6NcavEhGJKnU6sjfGTDfGFBlj1tR4fIQxJtsYk2uMeQDAWrvRWjuuxo+4AugEVAD5DbHwmmZlFTBx5moKin1YoKDYx8SZq5mVVdAYv05EJKrU9TTOq8CI6g8YY5zAC8BIIB243hiTfpT3pwGLrLW/AO46uaUe25S52fj8Aa5d4KLjbgOAryLAlLnZjfHrRESiSp2GvbV2IbC3xsNDgNzQkbwf+DtVR/C1yQf2hb4O1PYCY8wdxpilxpilu3btqsuyvqOw2Efb/YYfrHTx2KteRn3twhmoelxEJN7V5wJtJ2Bbte/zgU7GmDbGmBeBAcaYiaHnZgLDjTF/BBbW9sOstVOttRnW2ox27Wr9tO8xpaZ42Z1smTi+lOU9A/z4X24mz/CQcch7wj9LRCTW1OcCranlMWut3QNMqPFgKVDzPH6Dun94GhNnruZAswB/uaKcb9Ir+clnifzPnw15SXl0faQrTq+zMZcgIhKx6nNknw90qfZ9Z6Cwfss5eVcO6MRTV59JpxQvBiga7MbM6U7qbR3Z9uw2lvZbSvHC4nAtT0QkrOpzZJ8J9DTGdAMKgDHADQ2yqpN05YBO/x21PB/aj2lP9h3ZrPjBClLvSmXd9R6mfJWriKaIxI26Ri/fAhYBacaYfGPMOGttJXAPMBdYD7xjrV3beEs9ea1+2IrBqwbT+RedKfhrIcHL8mizzK+IpojEDWOtDfca/ktGRoZtrIrj0f8zj8vfgU57HHydXsnffljOwSTolOLlqwcubJTfKSLSFIwxy6y1GbU9F1F99saYUcCoHj16NNrvWNLSx7KxcNkiF5d94+KMzUm8cVE5mb0V0RSR2BVR3TjW2tnW2juSk5Mb7XekpnipTIBZ36/gkbE+drcM8j8fevh/s5MoLyxvtN8rIhJOETXsm8L9w9PwuqoimPntLI/fXMZ7F1XQe5ODJelLKHy5kEg8tSUiUh8RdRqnKRxO3RwuTOvY2stFz6ZxdotW5NyeQ87tORS9VUTaS2nMLdmrYjURiQlxd4H2WGzQsv3l7eTdn0elP8h73/fzcX8/NvTvH6/LyVNXn6mBLyIR6VgXaOPuNM6xGIch9Y5UBq8dTHbXINd+7uLXb3jotEvFaiIS3TTsa+Hp7OGZy0v5y6gy2hU7ePRVL1f8W8VqIhK94u6cfV2ltvKyON3Huq6l3PBFIld95SYjO4HZ14V7ZSIiJy6ijuyNMaOMMVNLSkrCvZQjqZ0DSfDXUeX8/sdlNC833PWig9xf5hIorbWpWUQkIkXUsG+KnH1d1SxW25XhxnzUndTbO5L/23wyz8pk34J9x/05IiKRQGmck7Bv/j6yb8+mLK+Mjnd0ZP2NXqZ8rWI1EQkvpXEaWKsLqorVuvyyC4Uvb8delke7pSpWE5HIpWF/kpxJTrpP6c6LE4IcTLTc+w8Pd36YSItSRTRFJPIojVNPS1r4WDq2qlht1KKqYrU3LypncR9FNEUkcujIvp5SU7wEnPDBeRVMGuujKCXIhNkefvVhEmX5ZeFenogIoGFfb9WL1QrbWR6/qYx3L64gbbODzPRMCv9aiA1G3kVwEYkvGvb1VDOimdray8XPpDN07RBaDG5BzoQcVv5wJaW5peFeqojEsYiKXla7ecntGzZsCPdy6s1ay47pO8i9Lxdbbjk0oTVPti0i/4AimiLS8I4VvYyoYX9YpOfsT1R5YTkLb1yFa8EhNp4SYPrIcvLbW7VoikiDUs4+zBJTE3nkkgO8cHkZbfc7eGSGlyv/5aLCp4imiDQNRS+bSGGJj4I+hIrV3Fz5tZvB2QlMH6lbIYpI49ORfRNJTfECcMgLL13m57fXlOHxw0Nvesj9RS6BQypWE5HGo2HfRKpHNAFWdw/w+AQ/FdekkP+7fDLPzGTfFypWE5HGodM4TaTmvW8Pp3GGD+hE8cJissdns/KilXQc35Fvb/byrIrVRKQBKY0TIQK+AJsf3czW/9tGidcy45JysnpWndpRakdE6kJpnCjg9Drp/nR3/nxXkP1JQX4208NdHyTS4pCK1USk/nQaJ8IsbeZj+S1w6WIXl3/tom+oWO2bdBWricjJi6gj+0i6LWG4HC5Wm/29qmK17a2D3PmRhwdmJVG2VcVqInJyImrYR9JtCcOlempne1vLkzeW8fbwCnptc5DZN5OCPxeoWE1ETlhEDXupvVht+FPpDF03hJbntGTD3RtYMWwFpTkqVhORulMaJ4pYa9kxYwd5P88j4AvQ7dFuLL/QwZQvchTTFJFjpnF0gTaKGGPoOLYjrYe3ZsPdG9j4wEb2nxLEOaIc2+E/978FNPBF5Dt0GicKJXZM5IyZZ/C3G4MkH4DJMzz8eKELV6VimiJSOx3ZR7HPO/v4ahxcP8/NqEVuBoWK1fJQTFNEvktH9lEsNcXLIS+8/CM//3dtGe5KePBND7cvTKLyYGW4lyciEUTDPopVj2muOT3Ar2/z8WVGgHO+MWSekcnez/aGeYUiEil0GieK1SxXa9PBy4CpaQw41Izs8dmsGr6KU8aeQvfnujNnc9F/lbDpIq5I/FD0MkYFygJseWwLW5/dSrCVk5eG+VjUveLI8ypXE4k9KkKLQ06Pk9OfPJ1BmYPY4arkzvfc3P1+IskHDaDUjki8iahhr26chtdiQAsm3VjKu+f76Z/n5IlpXs5dnQAWCouV2hGJFxE17NWN0zg6tPEy55wKJt3qo7BNkNs/TuS+dxJJt95wL01EmkhEDXtpHIdTO9vbWJ66sYzXLy6nZ6GTX/7RSf4f81WsJhIHNOzjQPVyNQx8e6GL8pndaH1+Crk/zSXr/CwOfXso3MsUkUakNE4cs9ay8/Wd5N6bS+BQgK6Tu5J1kZMp81SsJhKNlMaRWhljOOWWUxiyfghtr2jLpoc2ceDyHBK+Lcfyn2K1WVkF4V6qiNSThr3g7uCm7zt9efOmIC0OwqTXPFzzpQtXhSKaIrFCn6CVI/7ZycdX42HMPDeXfeNmUE5VsVquitVEop6O7OWI1BQvpR6YfqmfZ0f7SAjAQ296ufPLJCoPqFhNJJpp2MsR1YvV1nUN8uvbfMwbUsnZi6uK1fZ8uifMKxSRk6XTOHJEzWK1tu29DHwxjYFlzckel83qkavpcEsHejzXgzlbVawmEk0UvZQ6CZYH2fL4FrY+vZVASwcvX+Dj6+4VUFW1o2I1kQig6KXUmyPRQbffdGPQ0kFs91Zyxz/c3DNLxWoi0ULDXk5I837NmXRDKW8P83PWRidPvuzlvFUqVhOJdDpnLyfslNZePjnbx/Keldz6aSLjP0lk6LoEPr3OhHtpInIUEXVkr4rj6HA4tbOzteWZ68uYcUk53bc7uO+PTvL/kI8NRN51IJF4F1HDXhXH0aFmsVrOBS78s7rR5oJW5N6bS9Z5WRxap2I1kUiiNI40GGstRX8rYsPPNhA4EOC0X5/GiuEJKlYTaSJK40iTMMbQ4cYODFk3hHZXt2PzpM0cHJWDa72K1UTCTcNeGpy7vZv0t9J5/eYgzUph0useRs934VaxmkjYKI0jjWZeqo9F42D0fDcjl7gZsCGBV0aUk6NiNZEmpyN7aTSHi9VeGennmTE+HBYmvuVlwoIkKverWE2kKWnYS6OpXqy2/rSqYrV/nl3JkExDZt9M9sxRsZpIU9Gwl0ZTPaJpgHbtvGT8pQ+Dvh6IM9nJ6stWs+6mdfh3+8O9VJGYp+ilhEXQH2TLk1vY+uRWEpITKL6vDU/aQgpLFNEUOVmKXkrEcbgddHukG4OWDcLXwUGziTv48fQgyQeMIpoijUDDXsKq+ZnNeeQGH29dUE7fzVXFaj9YkYDPr4imSENS9FLCrmC/j/whsLxngFs/TeTWuYmcvT6BV0eUhXtpIjFDR/YSdqkpXgB2tbI8O6aMV4aX03WHg8df8bLtuW0qVhNpABr2EnbVI5oY+LJ/JY9NKMec05y8+/JY/r3lHFxzMLyLFIlyOo0jYVfz3repKV7uH53GRf1TKXq7iNz/zWXZwGWc9tBprBipYjWRk6HopUQ8/24/uT/LpehvRRS0C/LyiHI2pQYB3ftWpDpFLyWqudu6SX8zndduCeL1wcNveBgzz61iNZEToNM4EjXmd/Tx9Xi4boGbEZkuBm5wMn1EOdkqVhM5rog6stdtCeVYUlO8lCXCa8P9PHW9j6CBB/7u5a75SVSWqFhN5FgiatjrtoRyLNVTO9mnBpl0q4/Pzqlk8FLDkvQl7J69O8wrFIlcETXsRY6ltmK1IS/0YdDigbjauFhz+RrW3bAO/y4Vq4nUpDSOxISgP8jWZ7ay5TdbcLZ00vP5nizqXcGUzxTTlPihNI7EPIfbQdeHu5KRlYG3h5f1N65n2+hsfFvLdP9bETTsJcY069uMgV8NZM6PgvTa4uCJaV6GrUjAWMU0Jb4peikxxzgN753hY0Fnw9hPExk7N5Gh66ruf1uomKbEKR3ZS0xKTfGyK8UyZXQZ00eUc2qRg9+84uW6VV6ClcFwL0+kyWnYS0w6EtM0sLBfJQ+O87Hu9CAjP3GQdU4WB1epWE3ii4a9xKSaMc1mXTx0fSuN9LfTKdtSxrJBy9g0eRPB8iCzsgo49+l5dHtgDuc+PU8XcSUmKXopcadiTwW5P89l5+s7CZzu5rfnH2Bdh/98AlflahKtFL0UqcbVxkWf1/pw5sdncmBPOb981c31X7hxhz6LpdSOxCINe4lbbUa2YeKtpcwfUMnwpS4en+6lz+aqvxKFxUrtSGzRsJe41rqDl9cv8fPkDT6CDvjV215u/cTN6YnecC9NpEFp2EtcO5zayekS5OFbfcw52895qxN46M8uds3aFe7liTQYDXuJa9VTO5Uu+PrKBEpfP5UWnTysvWota69bi3+nitUk+imNI1KLYEWQbc9uY/Njm3E2d9Lj9z34pm+litUkoimNI3KCHC4Hpz10GhkrMkhKS+LbW76l4DoVq0n00rAXOYZmfZox4F8DmD0qSPetVcVqP1yuYjWJPipCEzkO4zTMTPfxZaeqYrWbP09kyHoVq0l00ZG9SB2kpnjZnWz5v+vKePnScjrvqipWG7PCS7BCxWoS+TTsReqgerHav8+s5MHxPtb0DDJ8roPlZy/nQNaBcC9R5Jg07EXqoGaxWvPOHrr9rTd93+tLeWE5ywYvY+NDGwmUBcK9VJFaKXopUk8VeyvIuy+PHa/uwJvmpehXbXh651ZFNKXJKXop0ohcrV30fqU3Z809i0P7K2h22zaGvRfA7VdEUyKHhr1IA2l9SWt+c6efLwZV8sNlCTwxzUvfTU5FNCUiRNSwN8aMMsZMLSkpCfdSRE7KZp+PNy/y8+SNZfgT4P53PIyb46ZkuyKaEl4RNeyttbOttXckJyeHeykiJyU1paotM7dzkMm3+vjwHD/fW5vAU9OT2PUPFatJ+ETUsBeJdkcimkBFAsw8v4Knx/nxdvaw9pq1rLlmDeU7ysO8SolH+gStSAM6nLqZMjf7SBrn7tFpXPiXjuT/Np9NkzdRPK+Y7s91Z/FZARWrSZNR9FKkCZVml5I9PpuSf5ewrluA6cPL2Z1c9XdQ976V+lL0UiRCJKUl0f/L/nx4eZBuBQ4en+blomUqVpPGp9M4Ik3MOAzv9/GxoJNh7KdubvpnImevT2C6itWkEenIXiQMUlO87G1pee7acqb+qJyOexw89qqX67NUrCaNQ8NeJAyqF6t9fUYlD44vZVWvIJd85mD5kOUcWK5iNWlYGvYiYVCzWK1FJy/d3+xN35l98e/ws2zIMjZO3EjAp2I1aRhK44hEmIp9FeTdn8eOaTvw9goVqxWpWE2OT2kckSjiauWi98u9Oevzszh4sIIW4/K54N0AieUqVpOTp2EvEqFaX1RVrDY3o4ILl1cVq525UcVqcnI07EUi2JZSH2/90M8TN5VR5ob73vUwXsVqchI07EUi2OFitbxOQSaP9fHB9/wMXZfA09OTKHq3iEi85iaRScNeJIJVL1arTID3v19VrObp4mHddetYe/VayrerWE2OT5+gFYlgtRWr3TM6jQv/3JH83+WzedJmlvRZQo/nerC4v4rV5OgUvRSJYqUbQsVqC0tY3y3AtEvK2Z2iYrV4peilSIxK6plE//n9+eDKIF0LHDwx3cvFSxMwQRWryXfpNI5IlDMOw6w0H192NPxkrpsbvwgVq41UsZr8h47sRWLA4WK1311Tzl8vK6PDPgePvurlxmVegn4Vq4mGvUhMqF6stqhvgAfHlbIiLchF/3SwbPAy9i/dH+4lSphp2IvEgJrFai07een5Rm/O+OAMKnZXsPzs5eT9vzwCpSpWi1dK44jEuMqSSvLuz2P7S9vx9vCS9nIaC1oe+k6cUzHN2KA0jkgcS0hOIG1qGv2+6IcNWlYMW8Gq8evZs9OHReVq8ULDXiROtLqwFYNXD+bf51nOy3LyxDQv/fKqPp2rmGbs07AXiSPOJCfTzi3l8ZvK8CXCz9/zcOfsRFqUQmGxYpqxTDl7kTiTmuJlIz4mj/Vx2SIXoxa56LspiY+vsFhrMcaEe4nSCHRkLxJnDsc0A0744LwKJo/1saeVZczfHay5Yg3lBSpWi0Ua9iJxpmZMk54eWn3Qi+6/7c6+f+5jSfoSCl8qxFrLrKwCzn16Ht0emMO5T8/TRdwopuiliBzhy/ORfXs2xfOLqcjw8tiQYra1qDzyvMrVIpuilyJSJ97uXvp90Y9eL/UisLqUX091M2JJVbEaKLUTzTTsReQ7jDGkjk9l4jgf67oGGDM/kYff8NBpV9WFW6V2opOGvYjUKqmLhz9cXc5fRpXRtqSqWO3Kf7vo0twb7qXJSdCwF5Fa3T88Da/byeL0qmK1Jb0DXPmVm8mveti/RMVq0UbDXkRqVT21cygJ5tzs4OAfUkkqNyw/Zzm59+WqWC2KKI0jIiekcn8lG3+1kcIXC/Gc7iHt5TS+TClVsVoEUBpHRBpMQssEev2lF/0X9Mc4DCsvXMnqcevZq2K1iKZhLyInJeUHKWSsymDh+ZZzVzh58mUv/TeoWC1SadiLyElzep28ck4pj91cxkEv3DvTw4QPE2lxSBHNSKMiNBGpl9QUL5vx8chPfFy62MXlX7vouzmJT0apWC2S6MheROqlerHa7O9VFavtam257h0Hq0etpmxbWbiXKGjYi0g91SxWMz08tJ7Vix6/70Hx/GIy+2ZS8GIBNhh5yb940uDRS2PM6cBDQLK19prQY8OA3wBrgb9baxcc62coeikSG3wbfWTfkU3xF8Ukn5/Mtl+k8Mz6TYpoNpJ6Ry+NMdONMUXGmDU1Hh9hjMk2xuQaYx4AsNZutNaOq/EjLHAQ8AD5J74JIhKNvKd76fd5P9KmpVG8/ACeazZz1txKTFARzaZW19M4rwIjqj9gjHECLwAjgXTgemNM+lHe/y9r7UjgV8CjJ7dUEYlGxhg63taRZ/63kjXdAoxe4Obh1z10KXIootmE6jTsrbULgb01Hh4C5IaO5P3A34ErjvL+UEEq+4DE2l5jjLnDGLPUGLN0165ddVq8iESP7KCP568q54XLy2i938HkGR6uXuiiaLcimk2hPhdoOwHbqn2fD3QyxrQxxrwIDDDGTAQwxlxtjPkr8Drwp9p+mLV2qrU2w1qb0a5du3osS0QiUWqKFwxk9gnw4PhSFvep5PJFbp54PYmSRSXhXl7Mq8+wry08a621e6y1E6y13a21T4UenGmtvdNaO/p4F2dFJDYdjmgCHPLCS5f5+dMYP20dLrLOzWLDvRsIHFKxWmOpz4eq8oEu1b7vDBTWbzkiEqsOp26qF6bdNDqN86d2YOPEjRT8oYA9H+yh10u9WNjGp2K1BlafYZ8J9DTGdAMKgDHADQ2yKhGJSVcO6FTr0O71p160H92e7HHZrLp4FWv7VbJvWDnW85/UzuH3y8mpa/TyLWARkGaMyTfGjLPWVgL3AHOB9cA71tq1jbdUEYllKd9PIWNlBl/+wDJ0lZMnpnkZmKNitYZSpyN7a+31R3n8Y+DjBl2RiMQtp9fJq0NLmdfVwW2fuPnp+x6WpFXyxsXlFKLUTn1EVBGaMWYUMKpHjx7hXoqIhElqipct+HjsljJGLnFxxVcu0reoWK2+Iqobx1o721p7R3JycriXIiJhUr1Y7aNzKph0q4+dbS3Xvutg9aWrKduqYrWTEVHDXkSkZrGao7uHNrN60eP5HhT/K1Ss9oKK1U6U7kErIlHDt9lHzh057Pt8H8nnhYrVslWsdpjuQSsiMcHb1ctZc88i7ZU09q08gPfazfT/VMVqdaFhLyJRxRhDx7FVxWoruwe49ks3k17zcOpOFasdi4a9iESlnICPP11Vzp+uLKPVQcPkGR5+rGK1o4qo6KWISF2lpngpKPaxNC3A+lN9jJnnZtQiN0NzXZRcVULyuUr1VRdRR/bGmFHGmKklJWrAE5Fjq1msNu1Hfp6/3k+bBBdZ389iw083UHmwMsyrjBxK44hI1JqVVfBfhWmX9ezApgc3UfCnAhK7JJI2NY2F7eOjWO1YaRwNexGJSSVflfDtuG/xZftYdFaAN4aVcchb9ZzX5eSpq8+MuYGv6KWIxJ3kc5PJWJHBgmGWIasdPDnNS0Z2/BaradiLSMxyepzMOLuUR39Sxr7mlntmebjn/USSDxoKi+MrtaM0jojEtNQUL1tDxWojlri46t8u+mzx8ull8VWspiN7EYlph1M7QQd8PLSCh2/1sb295cf/cLBq+Cp8m+PjCF/DXkRiWs1iNWd3D23f70XPF3qyf9F+Ms/IJP/5fGwg8sIqDSmi0jjV+uxv37BhQ7iXIyIxrmxrGTl35rD30720PKcladPS+LysOGpjmopeiogchbWWnW/sJPfeXCoOBPjwXD8fZvgJVAV3oiqmqeiliMhRGGM45eZTGLJuCKt7B7ligYvJMzyctqNqPMZKTFPDXkQEcHdw87tLS3n+qjJa+AyTXvNw7QIXrgpiIqap6KWISEhqipflvXx8e6qP0fPd/Gixm0E5Ccy+Ntwrqz8d2YuIhByOaZZ64JWRfp4Z7SPBwu1THeTcnUPl/ugtVtORvYhIyOGLsIfTOPv7J5Lw0550frec/D/ks2f2Hnr9tRdtRraptYQtki/iKo0jIlIHJd+UkD0um9J1pfh/1IIH03ax2xU48nwkpHaUxhERqafkoclkLM/gtIdPw/nJfia9mMjg9U4IHS9HempHw15EpI4ciQ66PdaNR37iY3fLIHd/6OGn7yeScqCqXyeSUzsa9iIiJyjYy8PjN5fx9jA/Z2xy8uQ0L+evTCA12RvupR1VRA173ZZQRKLB/cPTSEx08snZFfz6Nh9b2we57dNEHnm/Gb6NkXl0H1HD3lo721p7R3KybhQsIpGrernarlaWNyYYSh9sj+dbP5lnZrLt99sirlhNaRwRkQZSll9GzoQc9s7ZS4uzW9B7Wm8+9zddsZrSOCIiTcDT2cOZs8+kz5t98OX6WNI/k4U/X8eOPT4sUFDsY+LM1czKKmjytWnYi4g0IGMMHW7owJD1Q1iVbrn8SxePzPDSbXt4i9U07EVEGoG7nZvfjyzl91eX0dwHD7/u4br5LtxhKlZTXYKISCNJTfGyoqeP7FN9XDffzaVL3AzckMBH1zT9WnRkLyLSSA4Xq/kSYcYIP0+P8eEAxr/kIHtCNpUlTVespmEvItJIat7/9kC/RNwf9aTzfZ3Z/tJ2lvRdwp45e5pkLYpeioiEwf4l+8kel82hNYdof0N7NoxrxrNL8uoV0VT0UkQkwrQc0pJBywbR9ZGu7HynCMeojXT+2o+1jRPR1LAXEQkTh9tB18ldef7uIEXJQe6aXVWsZoINH9GMqDSOMWYUMKpHjx7hXoqISJNZ4fGx4ia4eFkCrQ4abOgwvCEjmhF1ZK9uHBGJR6kpXqwDPhtcydsXVHzn8YYSUcNeRCQeHY5oVud1Obl/eFqD/Y6IOo0jIhKPat77tjEK0zTsRUQiwJUDOjXq/Wt1GkdEJA5o2IuIxAENexGROKBhLyISBzTsRUTiQEQWoRljdgFb6vEj2gK7G2g54RQr2wHalkgVK9sSK9sB9duW06y17Wp7IiKHfX0ZY5YerfktmsTKdoC2JVLFyrbEynZA422LTuOIiMQBDXsRkTgQq8N+argX0EBiZTtA2xKpYmVbYmU7oJG2JSbP2YuIyHfF6pG9iIhUo2EvIhIHonbYG2NGGGOyjTG5xpgHanneGGOeDz2/yhgzMBzrrIs6bMswY0yJMWZF6M+kcKzzeIwx040xRcaYNUd5Ppr2yfG2JVr2SRdjzHxjzHpjzFpjzM9qeU1U7Jc6bku07BePMWaJMWZlaFsereU1DbtfrLVR9wdwAnnA6YAbWAmk13jNpcAngAGGAovDve56bMsw4KNwr7UO23I+MBBYc5Tno2Kf1HFbomWfdAQGhr5uAeRE8d+VumxLtOwXAzQPfe0CFgNDG3O/ROuR/RAg11q70VrrB/4OXFHjNVcAr9kq3wApxpiOTb3QOqjLtkQFa+1CYO8xXhIt+6Qu2xIVrLXbrbXLQ18fANYDNUvTo2K/1HFbokLov/XB0Leu0J+aaZkG3S/ROuw7AduqfZ/Pf+/0urwmEtR1neeE/sn3iTGmb9MsrcFFyz6pq6jaJ8aYrsAAqo4iq4u6/XKMbYEo2S/GGKcxZgVQBHxurW3U/RKtd6oytTxW8/+KdXlNJKjLOpdT1Xlx0BhzKTAL6NnYC2sE0bJP6iKq9okxpjnwD+Bea+3+mk/X8paI3S/H2Zao2S/W2gDQ3xiTArxvjDnDWlv9GlGD7pdoPbLPB7pU+74zUHgSr4kEx12ntXb/4X/yWWs/BlzGmLZNt8QGEy375LiiaZ8YY1xUDcc3rbUza3lJ1OyX421LNO2Xw6y1xcACYESNpxp0v0TrsM8Eehpjuhlj3MAY4MMar/kQuCV0RXsoUGKt3d7UC62D426LMeYUY4wJfT2Eqv22p8lXWn/Rsk+OK1r2SWiN04D11trnjvKyqNgvddmWKNov7UJH9BhjvMBFwLc1Xtag+yUqT+NYayuNMfcAc6lKs0y31q41xkwIPf8i8DFVV7NzgVLg1nCt91jquC3XAHcZYyoBHzDGhi7XRxJjzFtUpSHaGmPygclUXXiKqn0CddqWqNgnwLnAzcDq0PlhgAeBUyHq9ktdtiVa9ktHYIYxxknV/5DesdZ+1JgzTHUJIiJxIFpP44iIyAnQsBcRiQMa9iIicUDDXkQkDmjYi4g0MnOcYr2T+HmnGmM+C5XCrQt9oviYNOxFRBrfq/z3h6bq4zVgirW2D1X9WkXHe4OGvYhII6utWM8Y090Y86kxZpkx5l/GmN51+VnGmHQgwVr7eehnH7TWlh7vfRr2IiLhMRX4X2vtIOCXwJ/r+L5eQLExZqYxJssYMyX04axjispP0IqIRLNQmdv3gHdD7Q4AiaHnrgYeq+VtBdba4VTN7e9T1fq5FXgbGEtVlcRRadiLiDQ9B1Bsre1f84lQwVtthXWH5QNZ1tqNAMaYWVTd3OSYw16ncUREmliomnmTMeZaOHILwn51fHsm0MoY0y70/YXAuuO9ScNeRKSRhYr1FgFpxph8Y8w44EZgnDFmJbCWOt6hLtSD/0vgC2PMaqp671867hpUhCYiEvt0ZC8iEgc07EVE4oCGvYhIHNCwFxGJAxr2IiJxQMNeRCQOaNiLiMSB/w94QXYlk9+HnAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(t,realY(n0,t0,t),'m-')\n",
    "plt.scatter(t,w_RK4)\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "182c31bf",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
